Настройка конечно-элементной модели малого космического аппарата по результатам вибропрочностных испытаний

Повсеместное увеличение роли компьютерного моделирования при разработке изделий ракетно-космической техники, а также стремление к сокращению объёмов натурных вибропрочностных испытаний (ВИ) для решения задач прочности расчётным путём, с одной стороны, приводит к снижению стоимости этапа отработки прочности изделий, а с другой стороны, предъявляет повышенные требования к точности и детализированности конечноэлементных моделей (КЭМ).

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

Общее (укрупнённое) описание методики представлено в статье [1]. В ней рассматривается часть работы по созданию программного обеспечения (ПО) и методики автоматизированной валидации КЭМ по результатам ВИ, которая касается:

• •    особенностей выбора оптимальных параметров модели для настройки;

• •    реализации алгоритмов для анализа амплитудно-частотных характеристик (АЧХ) расчётных и испытательных откликов изделия.

Также в статье приведено сравнение результатов автоматизированной валидации КЭМ с результатами «ручной» настройки, представленными в статье [2].

В работе используется следующее ПО:

• •    Ansys Mechanical APDL для моделирования ВИ путём проведения гармонического анализа разработанной КЭМ изделия методом суперпозиции форм;

• •    утилита разработки отдела прочности РКК «Энергия» для эффективного (нештатного) получения результатов гармонического расчёта из Ansys ;

• •    DATADVANCE pSeven для организации процесса многопараметрической оптимизации в части математического аппарата, а также увязки всей системы в единую среду;

• •    программы на языке Python для реализации разрабатываемых в рамках работы алгоритмов в части сбора, подготовки и отображения данных, анализа АЧХ и вычисления значений целевой функции, которые используются в качестве критерия, по которому проводится оптимизация.

Следует отметить, что конкретное ПО может быть заменено на аналогичное при сохранении общей методики работ .

Описание изделия и его конечно-элементной модели

В качестве изделия, на примере которого в данной статье будут представлены результаты, выбран космический аппарат (КА) дистанционного зондирования Земли, разработанный (в части служебной платформы) и испытанный РКК «Энергия». В процессе создания КА были разработаны детальные динамические КЭМ аппарата и оснастки, с помощью которых были смоделированы натурные ВИ. После их завершения была проведена работа по валидации КЭМ, результаты которой представлены в статье [2].

Вышеупомянутый КА представляет собой конструкцию из несущих трёхслойных панелей с сотовым заполнителем, выполненных из АМг6. Приборное оборудование устанавливается на панелях через закладные элементы. Оснастка также выполнена из АМг6. Общая масса КА без оснастки — до 2 000 кг. Детальное описание конструкции представлено в статье [2].

Общий вид КЭМ КА с оснасткой и расположением контрольных точек, которые будут использованы в разделе статьи с примером работы методики, представлен на рис. 1. Модель состоит из 265 тысяч узлов. Диапазон модального расчёта составляет 0…200 Гц. В этом диапазоне расположены 670 собственных частот конструкции. Диапазон гармонического расчёта составляет 5…150 Гц.

В процессе работы, описанной в статье [2], подбор параметров производился «вручную». Было использовано 18 жёсткостных параметров, связанных с жёсткостью упругих связей (представленных в Ansys элементами типа COMBIN14 ), которые имитировали болтовые соединения в узлах конструкции, а также 22 диссипативных параметра, среди которых одно глобальное общеконструкционное демпфирование DMPRAT , два параметра, связанных с постоянным демпфированием в модели материалов DAMP и 19 параметров, связанных с частотно-независимым демпфированием в модели материалов DMPRAT .

Целью этапа работы, описанного в данной статье, была апробация алгоритмов и методики автоматизированной валидации КЭМ и улучшение результата, полученного «ручным» способом. Значения варьируемых параметров, достигнутые в статье [2], принимались как начальное приближение.

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

В процессе исследований были выявлены некоторые моменты, связанные с точностью настройки моделей:

• •    из-за ограничений КЭМ, связанных с невозможностью в компьютерной модели полностью воссоздать натурную схему жёсткостных и диссипативных связей, иногда требуется особый подход к моделированию стыков между узлами конструкции для достижения необходимой степени совпадения расчётных и испытательных откликов. Детальное описание этих подходов выходит за рамки данной статьи;

• •    диссипативные параметры, которые в рамках конечно-элементного анализа представляют собой реализованную тем или иным образом математическую абстракцию [3, с. 32], по необходимости могут принимать не только положительные, но и отрицательные значения, тем самым локально снижая уровень демпфирования элементов конструкции, заданный на более общем уровне.

Немаловажную роль в процессе автоматизированной валидации играет время расчёта одной итерации. Общее количество итераций в среднем может быть определено эмпирически, из расчёта, что на каждый варьируемый параметр приходится порядка 15–20 итераций. Таким образом, для рассмотренных выше 40 параметров в разработанной КЭМ общее количество итераций в процессе оптимизационного цикла составит приблизительно 600–800. Полная длительность гармонического расчёта методом суперпозиции форм обусловливается длительностью модального расчёта и зависит от исследуемого диапазона частот. Для решения задач динамической прочности относительно элементов крепления навесного оборудования диапазон частот составляет 0…~200 Гц. Так как методика автоматизированной валидации в данный момент требует настройки в части повышения точности решения, использовать крупные космические модули ввиду общей длительности расчётов неэффективно с точки зрения временных затрат, поэтому в данной статье рассматривается малый КА.

Благодаря тому, что (с некоторыми допущениями, описанными в работе [4, с. 32]) диссипативные параметры влияют на амплитуды резонансов КЭМ, а жёст-костные — в основном на частоты, становится возможным проводить валидацию по жёст-костным и по диссипативным параметрам раздельно, что:

• •    сокращает длительность решения (общее количество итераций) за счёт уменьшения количества варьируемых параметров в каждой из задач;

• •    для расчёта по диссипативным параметрам позволяет не проводить модальный анализ на каждой итерации, существенно сокращая её время.

В данной работе настройка КЭМ проводилась по обоим видам параметров (при этом методика позволяет проводить настройку как по жёсткостным, так и по диссипативным параметрам отдельно). В процессе натурных ВИ было использовано 65 трёхкомпонентных вибропреобразователей (ВиП); соответствующие им по расположению и ориентации контрольные точки были определены в КЭМ.

Алгоритмы автоматизированной валидации

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

Рис. 2. Укрупнённая схема вычисления общего коэффициента несоответствия TF: АЧХ — амплитудно-частотная характеристика; ВИ — вибропрочностные испытания; N_s — общее количество контрольных точек (датчиков); N_m — общее количество методов сравнения; M_f-1…M_f-Ns — общие коэффициенты несоответствия по всем методам для датчиков 1…N_s (рисунок создан автором)

Качество вычисления коэффициента несоответствия напрямую влияет на качество настройки КЭМ в процессе

Алгоритм сравнения откликов амплитудно-частотной характеристики по пикам оптимизации, так как:

• •    недостаточно точный анализ откликов АЧХ не позволяет определять различия между расчётными и испытательными АЧХ, что приводит к потере чувствительности от изменения варьируемых параметров КЭМ;

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

реакции на неважные, с точки зрения сравнения АЧХ, изменения или, наоборот, слишком малой реакции на важные изменения АЧХ.

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

Типовой вид графиков АЧХ (представляющих собой отклики вибрационного датчика) изображён на рис. 3. Здесь и далее по оси абсцисс показана частота в герцах, а по оси ординат — безразмерный коэффициент усиления (КУ). Он представляет собой отношение виброускорения конструкции в контрольной точке в единицах g ко входному виброускорению (в единицах g ), заданному с испытательного стенда или в расчёте. Таким образом, КУ показывает, во сколько раз входное вибрационное воздействие усиливается в определённой точке конструкции. Здесь и далее все графики АЧХ нормированы к единичному (1 g ) входному воздействию.

Рис. 3. Пример типовых амплитудно-частотных характеристик (АЧХ): a, б, в — отклики по осям X, Y, Z соответственно; =% — расчётная АЧХ; =% — экспериментальная АЧХ (рисунок создан автором)

Из приведённых графиков АЧХ видно, что расчётная модель не соответствует по своим откликам (по крайней мере, в указанной контрольной точке) испытательному изделию как по частотам резонансных пиков (РП), так и по их амплитудам. Задачей этапа сравнения откликов АЧХ является вычисление таких коэффициентов несоответствия по датчику M_{f -1}… M_f-Ns , которые бы выразили степень несоответствия графиков АЧХ одним числом — от нуля (идеальное совпадение кривых) до некоторого положительного значения (тем большего, чем хуже соответствие кривых).

Метод сравнения АЧХ по РП выделяет на графиках АЧХ зоны с резонансными пиками, вычисляет их характеристики, после чего ищет наилучшее соответствие между расчётными и испытательными резонансными пиками, которое в итоге послужит для вычисления коэффициента несовпадения по методу.

Схема вычисления коэффициента несовпадения по методу поиска РП представлена на рис. 4.

Прежде всего, отдельным алгоритмом, который будет рассмотрен ниже, выделяются РП на откликах. Затем создаётся матрица с числом столбцов по количеству найденных пиков на расчётном отклике i и числом строк по количеству найденных пиков на испытательном отклике j, которая заполняется коэффициентами TFij, характеризующими степень соответствия между собой каждой пары РП. TFij для пары пиков не вычисляется (элемент матрицы остаётся пустым), если абсолютное расстояние между двумя пиками превышает определённый глобальными параметрами порог Ffaraway (например, 10 Гц).

Элемент матрицы TF_ij определяется по следующим формулам:

при | f_i     j      faraway

TF_ij = NaN ;

при | f_i     j     faraway

^TF ij = ^ffreq ⁽ i’ j ) • ^k 1 ^{+ f} amp ⁽ i’ j ) ' ^k 2 ⁺

+ fabs-amp(i’ j) 'k3 + fnum(i’ j) ' k4’ где fi, fj — центральная частота расчётного и испытательного отклика соответственно; k1, k2, k3, k4 — весовые коэффициенты слагаемых как внешние параметры; ffreq(i, j) — абсолютная разница между центральными частотами РП, нормированная к порогу Ffaraway; famp(i, j) — абсолютная разница амплитуд РП, нормированная к максимальной амплитуде из них; fabs-amp(i, j) — абсолютная разница амплитуд РП, нормированная к максимальной амплитуде по датчику (во всём диапазоне частот); fnum(i, j) — абсолютная разница порядковых номеров РП, нормированная к максимальному числу пиков из двух откликов.

В табл. 1–3 представлены матрицы соответствия РП для датчика, отклики по осям X , Y и Z которого были представлены на рис. 3. Здесь зелёным цветом указаны наилучшие соответствия РП между собой, а серым — неоптимальные варианты.

Рис. 4. Схема вычисления коэффициента несоответствия по методу поиска резонансных пиков: АЧХ — амплитудно-частотная характеристика; A _расч[ x, y, z ] , A _исп[ x, y, z ] — расчётная и испытательная амплитудно-частотные характеристики соответственно; N _{пик.расч} , N _{пик.исп} — количество резонансных пиков (РП) на расчётной и испытательной АЧХ соответственно; PICK _расч , PICK _исп — наборы информации о расчётном или испытательном РП; M_pp — матрица соответствия расчётных и испытательных РП для текущих А _расч[ x, y, z ] и А _исп[ x, y, z ] (рисунок создан автором)

совпадений   по пикам

(размеры векторов) по откликам по осям X , Y и Z соответственно; W ( f ) — коэффициент веса пика, внешний параметр, зависящий от центральной частоты пика, выбранный таким образом,

Таблица 1

Матрица соответствия резонансных пиков, отклик по оси X

	РПрасч 1	РПрасч 2	РПрасч 3
РПисп 1	0,3842	—	—
РПисп 2	—	—	—

Примечание . РП_расч и РП_исп — расчётные и испытательные резонансные пики соответственно.

^Та^блица ² чтобы РП на низких

Матрица соответствия резонансных пиков, отклик по оси Y             частотах имели больший

	РПрасч 1	РПрасч 2	РПрасч 3	РПрасч 4	РПрасч 5	РПрасч 6	РПрасч 7	РПрасч 8
РПисп 1	—	—	—	—	—	—	—	—
РПисп 2	—	—	—	—	—	—	—	2,3137

вес, нежели на высоких, для повышения приоритета в диапазоне частотного отклика элемен-

	тов крепления навесного
Матрица соответствия резонансных пиков, отклик по оси Z	Таблица 3 оборудования (приблизительно 20…150 Гц).

	РПрасч 1	РПрасч 2	РПрасч 3	РПрасч 4	РПрасч 5	РПрасч 6	РПрасч 7
РПисп 1	1,9396	1,1664	—	—	—	—	—
РПисп 2	—	—	—	1,9310	1,3975	—	—
РПисп 3	—	—	—	—	2,2042	—	—
РПисп 4	—	—	—	—	—	2,5500	—

Для откликов по осям X и Y при текущих значениях варьируемых параметров КЭМ соответствие можно установить лишь по одной паре пиков. Для отклика по оси Z и ситуаций, когда существует несколько возможных соответствий, алгоритм поиска лучших совпадений следующий:

• •    выбирается наименьшее значение коэффициента TF_ij , в данном случае tf = 1,3975, между расчётным пиком РП_расч 5 и испытательным пиком РП_исп 2;

• •    принимается, что РП_расч 5 соответствует РП_исп 2 и остальные элементы матрицы с их участием обнуляются;

• •    аналогичным образом находится соответствие между РП_расч 2 и РП_исп 1, а также РП_расч 6 и РП_исп 4.

Таким образом, для отклика по оси Z определено совпадение трёх РП.

Общий коэффициент FM_pp определяется по следующей формуле:

W < f )       W < f )       W ( f )

{ tf } '        ⁺ { tf } '        ⁺ { f } '

x N y N z N tfx                      tfy                     tfz

FM pp

В части выделения РП из АЧХ применяется алгоритм, схематично изображённый на рис. 5.

Для поиска РП кривая АЧХ упрощается путём поиска среднего значения на участке заданной во внешнем параметре ширины, тем самым исключается влияние мелкого шума по амплитуде. Затем выявляются зоны, в которых присутствуют положительные экстремумы (пики на графике).

Рис. 5. Схема выделения расчётных пиков (РП) на амплитудно-частотной характеристике (АЧХ): A _расч[ x, y, z ] , A _исп[ x, y, z ] — расчётная и испытательная амплитудночастотные характеристики соответственно; a_thr — средняя амплитуда по отклику, ниже которой РП отбрасываются; a_mid — средняя арифметическая амплитуда по отклику; F_CUTOFF — настраиваемый весовой коэффициент (рисунок создан автором)

После чего анализируется оригинальная АЧХ в зонах, найденных на предыдущем шаге, и для каждой зоны по частоте с пиком на упрощённом графике выбирается максимальная точка на оригинальном графике, которая будет представлять собой вершину резонансного пика. Упрощённые кривые на рис. 3 представлены пунктирными линиями, а сами центральные частоты резонансных пиков — вертикальными линиями. Полупрозрачные зоны отображают условную ширину РП.

Алгоритмы прямого сравнения амплитудно-частотной характеристики и вычисление общего коэффициента несоответствия

Метод прямого сравнения АЧХ вычисляет коэффициент несовпадения интегрально, как сумму разницы средних значений АЧХ на всём заданном частотном диапазоне с некоторым шагом.

Схема вычисления коэффициента несовпадения по прямому методу представлена на рис. 6.

Суть метода состоит в сканировании частотного диапазона с переменным шагом по частоте: чем выше частота, тем больше шаг, для исключения влияния отдельных пиков на высоких частотах, которые не интересны с точки зрения прочности навесного оборудования [1, с. 47], [5, с. 74]. Абсолютная разница между средними значениями амплитуд на шаге прежде всего нормируется по максимальной амплитуде на текущем участке для исключения влияния на величину коэффициента абсолютного значения амплитуды. Общая сумма представляет собой коэффициент несовпадения откликов по методу, который нормируется по частоте для исключения влияния различных по длине диапазонов и по максимальной амплитуде откликов для сглаживания неравномерности между максимальными амплитудами откликов одного датчика.

Рис. 6. Схема вычисления коэффициента несоответствия по методу прямого сравнения: АЧХ — амплитудно-частотная характеристика; КУ — коэффициент усиления; F _НД , F _КД — частота начала и конца диапазона соответственно; f, f_s — текущая частота и шаг по частоте соответственно; m _расч , m _исп — средняя величина КУ на исследуемом участке (шаге); MA _расч , MA _исп — максимальная величина КУ во всём диапазоне текущего отклика; d, D — величина разницы между расчётными и испытательными средними КУ на участке и их нормированная разница соответственно; S — накопленная сумма нормированной разницы величины средних испытальных и расчётных КУ на участке; FM _dir[ x, y, z ] — коэффициент несовпадения текущих расчётных и испытательных откликов; FM _dir – общий коэффициент несовпадения по всем трём откликам для метода DIRECT (рисунок создан автором)

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

После того как для всех участвующих в анализе датчиков (контрольных точек) найдены коэффициенты несоответствия по всем используемым методам (прямое сравнение и поиск пиков в данном случае), происходит их объединение в общий коэффициент по датчику по следующей формуле:

WMf = W• (FMdir • kdir + FMpp • kpp), где W — пользовательский весовой коэффициент конкретного датчика, который выбирается исходя из экспертной оценки вклада датчика в общий результат; kdir, kpp — весовые коэффициенты соответствующих методов для исключения разбалансировки алгоритма в сторону одного из методов определения несоответствия.

Сумма WMf по всем датчикам, поделённая на количество датчиков, является искомым общим коэффициентом несоответствия, который передаётся в pSeven в качестве результата одной итерации расчёта и подлежит минимизации путём подбора варьируемых параметров.

Пример работы алгоритмов автоматизированной валидации для космического аппарата дистанционного зондирования Земли

В качестве примера работы системы автоматизированной валидации рассмотрим показанный ранее КА ДЗЗ. В качестве варьируемых параметров использовались:

• •  18 диссипативных параметров

материалов типа DMPR ;

• •    один диссипативный параметр материала типа DAMP ;

• •    один общеконструкционный диссипативный параметр типа DMPRAT ;

• •  10 диссипативных параметров

CV1 для элементов типа COMBIN14 ;

• •    восемь жёсткостных параметров (линейных жёсткостей «пружин») для элементов типа COMBIN14 .

Прочие жёсткостные параметры (кроме выбранных восьми) не варьировались, а значения их были приняты в соответствии с лучшими значениями, подобранными «вручную» с результатами, представленными в статье [2]. Полное исключение из оптимизации жёсткостных параметров позволило бы сильно сократить общее время расчётов благодаря отсутствию необходимости при каждой итерации проводить модальный анализ. Вместо этого в Ansys Mechanical APDL возможно применение гораздо более быстрого пересчёта элементных матриц, основываясь на единожды полученных результатах модального анализа (собственных форм и частот). Анализ проводился по 65 контрольным точкам (трёхкомпонентным датчикам). Для достижения приведённых результатов цикл оптимизации составил ~1 000 итераций с общим временем решения 140 ч. При этом 120 ч было потрачено на непосредственную работу ПО Ansys (расчётная станция позволяла без потери производительности запускать задачу в два параллельных потока по 16 ядер с тактовой частотой 2,6 ГГц), так что без этого общее линейное время составило бы 240 ч.

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

Исходная модель была доработана в части замены некоторых, выбранных экспертной оценкой, связей типа COUPLING элементами типа COMBINE14 для возможности задания в них диссипативных свойств CV1 . Для доработки были выбраны элементы конструкции, влияющие на отклики в середине исследуемого частотного диапазона, приблизительно 40…80 Гц. Для сохранения общей картины распределения жёсткостей в дорабатываемых узлах «бесконечная» жёсткость связей COUPLING была заменена на высокую, но конечную жёсткость в COMBINE14 .

На рис. 7–9 приведены графики АЧХ для некоторых наиболее показательных датчиков до и после автоматизированной настройки КЭМ. Расположение контрольных точек показано на рис. 1. (см. с. 7)

5    15   25   35   45   55   65   75   85   95 105 115 125 135 145

Частота, Гц а)

5    15   25   35   45   55   65   75   85   95 105 115 125 135 145

Частота, Гц б)

Частота, И'щ ^

Рис. 7. Амплитудно-частотные характеристики датчика № 6: а, б и в — отклики по осям X, Y и Z соответственно; ^— — испытания; ^— — до настройки конечно-элементной модели (КЭМ), TF = 1,25; ^— — после настройки КЭМ, TF = 1,11; (рисунок создан автором)

Рис. 8. Амплитудно-частотные характеристики датчика № 35: а, б и в — отклики по осям X, Y и Z соответственно; ^— — испытания; ^— — до настройки конечно-элементной модели (КЭМ), TF = 1,16; ^— — после настройки КЭМ, TF = 1,00 (рисунок создан автором)

Рис. 9. Амплитудно-частотные характеристики датчика № 46: а, б и в — отклики по осям X, Y и Z соответственно; ^— — испытания; ^— — до настройки конечно-элементной модели (КЭМ), TF = 1,61; ^— — после настройки КЭМ, TF = 1,14 (рисунок создан автором)

Среднее снижение коэффициента несоответствия по всем используемым в расчёте датчикам составило приблизительно 20%. При этом можно отметить, что разницы частот и амплитуд первых собственных резонансов (на частотах до 20 Гц) составляют менее 20%, что можно считать приемлемым значением расхождения для первых собственных тонов [5, с. 79], [6, с. 24] и [3, с. 35].

Также в «среднем» диапазоне 30…100 Гц можно наблюдать более выраженные РП с тенденцией к смещению в сторону соответствующих испытательных РП. На частотах свыше 100 Гц наблюдается плохое совпадение откликов, что может говорить об отсутствии среди варьируемых параметров тех, что связаны непосредственно с высокочастотными откликами конструкции или недостаточной подробностью КЭМ.

Выводы

Разрабатываемая методика показала работоспособность на относительно малых по размеру КЭМ (до 500 тысяч узлов).

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

Для решения задачи настройки крупных КЭМ рекомендуется разработка более эффективных подходов, которые позволят сократить общее расчётное время.

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

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