Full Text

Введение Конечно-элементные модели (КЭ-модели) широко используются для проведения статистических и динамических расчетов во многих областях техники. Однако такие модели в ряде случаев содержат неизбежные погрешности моделирования, обусловленные дискретизацией модели, неточностью задания свойств материалов, геометрических характеристик и граничных условий [1]. Для устранения погрешностей моделирования были разработаны различные методы коррекции КЭ-моделей применительно к разным конструкциям: топливным бакам ракет-носителей [2], соплу ракетного двигателя [3], композитному самолету [4], архитектурным сооружениям [5], пластинчатым теплообменникам [6]. В основе этих методов лежит минимизация разницы между ключевыми характеристиками реальной конструкции и параметрами КЭ-модели. В качестве таких параметров могут выступать, например, частоты собственных колебаний и отклик конструкции на динамическое воздействие [7]. Кроме того, методы коррекции нередко используются для оценки поврежденностей различных конструкций: дамб [8], мостов [9, 10], контейнеров для хранения радиоактивных отходов [11]. Обзор подходов к решению проблемы мониторинга поврежденности приведен в работе [12]. Следует отметить, что погрешность описания механических свойств реальной конструкции методом конечных элементов зачастую обусловлена особенностями моделирования узлов сопряжения, например болтовых и заклепочных соединений. Однако описание нелинейного поведения таких узлов посредством построения нелинейной КЭ-модели в ряде случаев может оказаться вычислительно затратным. Во избежание этого такие узлы сопряжения могут быть заменены специальными элементами с эквивалентными жесткостными и диссипативными характеристиками. Так, в работах [13, 14] болтовые узлы сопряжения двухбалочной модели были заменены гистерезисными элементами Iwan’a, параметры которых были определены методом коррекции. Известные методы коррекции могут быть разделены на две категории: стохастические и детерминированные. В основе стохастических методов лежит представление о том, что экспериментальные данные являются случайными и содержат неизбежные ошибки измерения, обусловленные как объективными факторами (температура, влажность, вибрации, оборудование), так и субъективными (опыт проведения экспериментов) [15]. В зависимости от типов ошибок измерения были разработаны различные методы коррекции: метод симуляции Монте-Карло [16, 17], методы пертурбаций [18, 19], методы случайных поверхностей отклика [20, 21] и Байесовcкие методы [22, 23]. Однако стохастические методы коррекции на несколько порядков вычислительно затратнее детерминированных методов. Заметим, что вопрос выбора параметров коррекции для стохастических методов не является однозначным [24]. Детерминированные методы коррекции обычно сводятся к итерационной процедуре коррекции параметров КЭ-модели. Среди итерационных методов наиболее популярным и интуитивно понятным является подход, основанный на создании матрицы чувствительности [25-28]. В основе этого метода лежит проблема минимизации целевой функции, равной сумме квадратов разностей между измеренными в эксперименте данными и соответствующими данными, полученными с помощью модели. Удобство такого подхода заключается в том, что в целевой функции одновременно могут фигурировать разнородные параметры, такие как собственные частоты и отклики конструкций на динамическое воздействие. Чтобы уравновесить вклад этих данных в результирующую функцию, авторы работы [29] вводят весовые коэффициенты, значения которых определяют методом анализа иерархий. Нередко применение методов коррекции приводит к тому, что результирующая система уравнений является плохо обусловленной [30]. Для борьбы с этой проблемой существуют две техники регуляризации, наиболее часто используемые исследователями: регуляризация Тихонова [31, 32] и сингулярное разложение [33]. Особенности применения обоих подходов на примере модели бетонного сооружения показаны в работе [34]. В данной работе предлагается детерминированный метод коррекции КЭ-моделей конструкций по результатам модальных испытаний. Матрица масс считается определенной точно и в ходе коррекции не изменяется. 1. Метод коррекции конечно-элементной модели Пусть задана конечно-элементная модель исследуемого объекта в виде матрицы жесткости и матрицы масс , тогда собственные числа (где - частота собственных колебаний) и формы колебаний определяются из решения обобщённой проблемы собственных значений: . (1) Для коррекции динамических свойств КЭ-модели будем вносить изменения только в матрицу жесткости путем добавления к исходной матрице матрицы жесткости корректирующей КЭ-модели . Для модели из объемных элементов предлагается использовать в качестве корректирующей КЭ-модели невесомую ферменную конструкцию, которая строится на узлах исходной КЭ-модели. Заметим, что жесткость корректирующих элементов может быть отрицательной для изменения упругости конструкции в сторону уменьшения. Жесткость каждого такого «ферменного» элемента является параметром коррекции, подлежащим определению из решения задачи безусловной минимизации целевой функции. На рис. 1 приведена принципиальная схема, иллюстрирующая физическую сторону предлагаемого метода на примере четырехузлового КЭ в виде пространственного тетраэдра. В данном случае добавлено 6 ферменных элементов (по сути, это 6 пружин) для коррекции упругости суммарной модели. Таким образом, корректирующая модель образует «каркасную» структуру, которая подходит для коррекции моделей из объемных КЭ-элементов. В случае, когда корректируемые динамические свойства модели существенно зависят от изгибных и крутильных жесткостей её балочных или оболочечных элементов, необходимо вводить корректирующую модель из балочных элементов. Рассмотрение этого случая выходит за рамки данной статьи. Рис. 1. Создание корректирующей КЭ-модели на примере четырехузлового КЭ в виде пространственного тетраэдра Fig. 1. Creating the updating finite element model illustrated by the example of a tetrahedrical finite element Матрица жесткости корректирующего элемента строится как ферменный конечный элемент и имеет следующий вид: , (2) где - направляющие косинусы; - неизвестная корректирующая жесткость; - число корректирующих элементов. Матрица (2) симметричная, поэтому не нарушает симметрию матрицы жесткости корректируемой модели при суммировании с ней. Для автоматического формирования набора корректирующих элементов используется портрет матрицы жесткости корректируемой конструкции. Так, если в матрице линейные степени свободы двух узлов «связаны» между собой внедиагональными элементами, то добавляется корректирующий элемент между этими узлами. В общем случае число корректирующих элементов определяется количеством связей между узлами в матрице, но оно может быть уменьшено посредством выбора областей коррекции, например элементов конструкции с наибольшей неопределенностью физических и геометрических характеристик. Кроме того, число независимых корректирующих жесткостей может быть уменьшено посредством введения зависимостей между элементами. Например, могут быть учтены геометрические особенности конструкции: наличие плоскостей симметрии и конструктивно идентичных элементов. Также можно проводить поэтапную коррекцию, в ходе которой в качестве целевых принимаются различные группы частот. В зависимости от геометрии модели и типов конечных элементов может получиться так, что при добавлении (2) к матрице жесткости исходной конструкции появятся новые ненулевые элементы (изменится портрет разреженной матрицы). Такая ситуация возникает, например, когда есть геометрическая симметрия относительно некоторого узла и в исходном портрете соответствующие элементы при ассемблировании матрицы жесткости в сумме дают ноль и не попадают в портрет. Поэтому перед началом процедуры коррекции портрет матрицы жесткости дополняется портретом корректирующей матрицы. Рассмотрим алгоритм метода коррекции на примере одного целевого значения для упрощения индексных обозначений, а затем обобщим формулы для случая коррекции по нескольким значениям собственных частот. Пусть в результате коррекции необходимо достигнуть целевого значения . Необходимо найти такую , чтобы достигнуть значения при решении следующей обобщенной проблемы: . (3) Так как в (3) собственный вектор зависит от матрицы , то использовать (3) напрямую как нелинейное матричное уравнение для поиска корректирующих жесткостей слишком затратно с вычислительной точки зрения: придется многократно решать обобщенную проблему. Поэтому для поиска корректирующей матрицы предлагается следующий итерационный алгоритм. Пусть решена обобщенная проблема для исходной модели , (4) то есть найдены и . Здесь и далее верхний индекс в круглых скобках обозначает номер итерации. На первой итерации ищем такое , чтобы удовлетворить следующему уравнению: . (5) Здесь необходимо отметить, что является функцией от корректирующих жесткостей . Число корректирующих элементов в общем случае больше размерности матрицы жесткости , а следовательно, больше числа уравнений в (5). Поэтому для решения (5) необходимо устремить к нулю следующую скалярную величину (невязку): , (6) где . В результате минимизации невязки (6) по корректирующим жесткостям получим и, решив новую обобщённую проблему , (7) найдем и . Это первая итерация метода. Аналогично (5) составим следующее уравнение для поиска : , (8) представим , тогда (8) перепишется в следующем виде: (9) где . Аналогично (6) устремим к нулю невязку . (10) Это была вторая итерация метода коррекции. Обобщим формулы на случай нескольких целевых собственных значений. Пусть их число равно . Таким образом, на каждой итерации метода коррекции минимизируется следующая целевая функция: (11) где - номер итерации; - весовой коэффициент i-го тона. Для минимизации (11) используется метод сопряженных градиентов. С вычислительной точки зрения важно найти аналитические выражения для вектора-градиента целевой функции. Для этого найдем производные от целевой функции по неизвестным корректирующим жесткостям (здесь и далее опустим номер итерации в круглых скобках): , . (12) Пусть - компоненты вектора собственной формы , тогда частные производные из (12) примут вид , , . (13) Для вычисления (13) нет необходимости «прокручивать» два цикла (суммы), так как ненулевыми компонентами в общем случае без учета симметрии являются максимум 36 элементов в (2), поэтому (13) можно заменить на циклы по элементам (2). Рассмотрим такую замену на примере одного корректирующего элемента. Для удобства изложения обозначим произведения направляющих косинусов как и введем локальную нумерацию элементов, входящих в произведения (13): (14) На основе выражений (14) можно вычислить производные . Рассмотрим другой вариант вычисления производных. Для этого воспользуемся тем фактом, что симметричные матрицы хранятся в виде верхнего или нижнего треугольника, причем в разреженном формате хранения матрицы её элементы отсортированы в порядке возрастания номеров строк и столбцов. Для определенности рассмотрим построчное хранение разреженной матрицы: . (15) В (15) введена другая нумерация произведений направляющих косинусов, и для удобства отрицательный знак внесен в значения соответствующих элементов. В общем виде один корректирующий конечный элемент с учетом симметрии имеет максимально 21 ненулевой элемент. В зависимости от конкретных значений направляющих косинусов ненулевых элементов в (15) может быть меньше. Для каждого корректирующего элемента предлагается хранить все произведения направляющих косинусов в виде последовательности , где - количество ненулевых элементов -го корректирующего конечного элемента. Кроме того, для каждого такого элемента хранится порядковый номер элемента в портрете разреженной матрицы жесткости, по которому можно определить строку , столбец (от 1 до ) и соответствующие номера степеней свободы , (принимают значения от 1 до 3, соответствуют перемещениям вдоль осей глобальной системы координат). Тогда формула для вычисления производной (13) примет вид , (16) где - символ Кронекера. Если сравнить вычислительные затраты на (14) и (16), первый имеет примерно в 2 раза меньше умножений, но для его реализации требуется дополнительная память на массив, содержащий для каждого корректирующего элемента. В процессе итераций при «движении» к целевым значениям частот возможно изменение порядка частот, особенно при наличии в спектре близких и кратных частот, поэтому после каждой итерации необходимо находить соответствие новых тонов колебаний старым. Обычно это можно сделать, используя критерий модального соответствия (MAC-критерий): , однако при наличии кратных частот колебаний MAC-критерий может дать близкие численные значения сразу для нескольких форм колебаний, что не позволяет установить однозначное соответствие между формами колебаний до и после коррекции. Для решения этой проблемы предлагается находить собственные частоты и формы колебаний на каждом шаге коррекции методом итераций в подпространстве [35], используя в качестве начальных приближений векторы, найденные на предыдущем шаге. Рассмотрим более общий алгоритм, который в частном случае вырождается в представленный выше. Представим , тогда (5) перепишется в следующем виде: . (17) Вычтем (4) из (17), составим целевую функцию для поиска : где . (18) Решив задачу минимизации (18), получим и решим новую обобщенную проблему: . (19) При этом , если и . Сформируем целевую функцию для поиска второго приближения . Невязка будет иметь следующий вид: ; . Умножим слева на транспонированный собственный вектор, получим (20) Решив задачу минимизации (20), будем иметь . Решим новую обобщенную проблему: . При этом , если и . Сформируем целевую функцию для поиска третьего приближения : . В итоге аналогично (20) получим (21) Таким образом, проанализировав (20) и (21), для поиска корректирующей матрицы жесткости имеем следующую целевую функцию: (22) Выберем следующую нормировку для собственных векторов : . После решения обобщенной проблемы нормируем - сначала , для установления соответствия тонов колебаний до и после коррекции по MАС-критерию, а затем , тогда , . Таким образом, целевая функция примет следующий вид: . (23) Выражения для производных от целевой функции, входящие в градиент по типу (16), имеют следующий вид: (24) где - компоненты вектора , а - компоненты вектора . Таким образом, в целевой функции (23) содержатся собственные векторы для двух обобщенных проблем: исходная и достигнутая в процессе итераций проблема для скорректированной матрицы для поиска следующего приближения корректирующей матрицы . Если обновлять «базовую точку» итерационного процесса на каждой итерации, то получим первый алгоритм с целевой функцией (10) при условии выбора нормировки собственных форм к единичным обобщенным массам: . Сходимость представленного итерационного алгоритма существенно зависит от геометрии и физического соответствия расчетной модели реальной конструкции, поэтому она определяется практическим путем в каждом конкретном случае на основе требований, предъявляемых к расчетной модели летательного аппарата. Так, например, если назначением модели является оценка эксплуатационной нагруженности летательного аппарата, то, как правило, коррекция расчетной модели призвана обеспечивать соответствие с заданной точностью расчетных значений собственных частот низших тонов колебаний их экспериментальным значениям. В то же время решение проблемы аэроупругой устойчивости самолетов требует сохранения в расчетной модели собственных частот таких тонов, которые оказывают первостепенное влияние на исследуемое явление. И степень соответствия модели реальной конструкции здесь определяется ограничениями максимальной скорости летательного аппарата. 2. Апробация метода Алгоритм предложенного метода коррекции реализован в виде программы на языке Fortran. Программа была успешно протестирована на простых КЭ-моделях, составленных из оболочечных и объемных элементов. Для апробации метода на примере реальной конструкции была выбрана динамически подобная модель самолета Ту-204 (ДПМ), выполненная по отсечно-балочной схеме (рис. 2). Габаритные размеры: размах крыла 3172 мм, длина фюзеляжа 3462 мм (масштаб моделирования 1:10). Масса с датчиками и кабелями 50,5 кг. Упругие характеристики фюзеляжа и крыла смоделированы лонжеронами, расположенными по оси жёсткости каждого агрегата. Топливо, коммерческая нагрузка, оборудование, носовая и основные опоры имитировались жёсткими сосредоточенными грузами. Для проведения экспериментального модального анализа ДПМ была вывешена на упругой подвеске малой жесткости. Рис. 2. Общий вид ДПМ на упругой подвеске Fig. 2. General view of the dynamically-scaled model with the flexible suspension Создание 3D-модели (рис. 3, a) производилось по натурной ДПМ. Геометрические модели крыльевых отсеков не создавались, так как они, в силу конструктивного исполнения ДПМ, являются точечными массами для упругого крыла - они не вносят жесткостей и не соприкасаются между собой при малых амплитудах колебаний. В силу последнего обстоятельства отсеки и хвостовое оперение при создании КЭ-модели также моделировались дискретными массами с инерционными характеристиками, определенными по соответствующим геометрическим моделям. Число степеней свободы в результирующей КЭ-модели (рис. 3, б) составило 751871. Коррекция КЭ-модели проводилась по шести наборам экспериментально определенных частот: 1) по 7 тону - симметричный изгиб крыла I тона (СИКр1); 2) с 7 по 8 тон - СИКр1, антисимметричный изгиб крыла I тона (АСИКр1); 3) с 7 по 9 тон - СИКр1, АСИКр1, горизонтальный изгиб фюзеляжа I тона (ГИФ1); 4) с 7 по 10 тон - СИКр1, АСИКр1, ГИФ1, вертикальный изгиб фюзеляжа I тона (ВИФ1); 5) с 7 по 10 и 12 тону - СИКр1, АСИКр1, ГИФ1, ВИФ1, симметричный изгиб крыла II тона (СИКр2); 6) с 7 по 10, 12 и 15 тонам - СИКр1, АСИКр1, ГИФ1, ВИФ1, СИКр2, вертикальный изгиб фюзеляжа II тона (ВИФ2). а б Рис. 3. Геометрическая (a) и конечно-элементная (б) модели ДПМ Fig. 3. Geometrical (a) and finite element (b) models of the dynamically-scaled model Итерационный процесс коррекции по предложенному алгоритму считался завершенным при достижении целевых значений частот с точностью 0,0001 %. Результаты коррекции сведены в таблице. Распределения изменений узловых жестокостей по всем линейным степеням свободы КЭ-модели до и после каждой из коррекций приведены на рис. 4, a-е. Синяя цветовая гамма на рисунках соответствует областям понижения исходной жесткости (корректирующая КЭ-модель имеет отрицательную жесткость в этих областях), а красная - областям повышения исходной жесткости. Черному цвету соответствуют области, жесткость которых в ходе коррекции осталась неизменной. Заметим, что коррекция по одному тону (набор 1) приводит к снижению жесткости практически во всех областях (рис. 4, a) с измененной жесткостью. Добавление в набор других тонов приводит к появлению областей с положительным изменением жесткости в узлах (рис. 4, б-е). При этом для каждого набора распределение корректирующих жесткостей уникально. Вследствие изменения упругости корректируемой модели происходит изменение частот и форм колебаний. Нормированные к массе формы колебаний до (серым цветом) и после коррекции (зеленым цветом) по шестому набору частот приведены на рис. 5. Из рисунков видно, что формы колебаний до и после коррекции хорошо коррелируют между собой, разница между формами колебаний по MAC-критерию до и после коррекции составила: 0,02, 1,12, 1,74, 0,87, 0,18 и 4,09 % соответственно. Результаты коррекции Results of the model updating Тон Частота, Гц Погрешность до и после коррекции, % эксперимент исходная модель до после 1 2 3 4 5 6 СИКр1 3.440 3.492 1.5 0.0 0.0 0.0 0.0 0.0 0.0 АСИКр1 4.873 4.955 1.7 -0.6 0.0 0.0 0.0 0.0 0.0 ГИФ1 5.440 5.728 5.3 4.7 4.7 0.0 0.0 0.0 0.0 ВИФ1 5.730 5.968 4.2 3.5 3.3 2.6 0.0 0.0 0.0 СИКр2 9.300 9.131 -1.8 -4.4 -3.2 -3.9 -4.4 0.0 0.0 ВИФ2 14.18 14.77 4.2 3.7 3.8 3.4 2.0 3.3 0.0 a б в г д е Рис. 4. Распределения изменений узловых жесткостей при коррекции: по наборам целевых частот от 1 до 6 (a-е) Fig. 4. Alterations of nodal stiffnesses during updating: by the sets of the target frequencies from 1 to 6 (a-f) a б в г д е Рис. 5. Формы колебаний до и после коррекции: a - СИКр1; б - АСИКр1; в - ГИФ1; г - ВИФ1; д - СИКр2; е - ВИФ2 Fig. 5. Modeshapes before and after updating: a - the first symmetric wing bending; b - the first asymmetric wing bending; c - the first horizontal fuselage bending; d - the first vertical fuselage bending; e - the second symmetric wing bending; f - the second vertical fuselage bending Заключение Предложен метод коррекции расчетных динамических моделей, который состоит в том, что к исходной матрице жесткости добавляется матрица жесткости корректирующей КЭ-модели, при этом матрица масс считается определенной точно и не изменяется в ходе коррекции. В рамках данной работы корректирующая модель представляет собой невесомую ферму, построенную на узлах исходной модели. Жесткости корректирующих элементов вычисляются итерационно из решения задачи минимизации целевой функции, которая является взвешенной суммой квадратов разности текущих и целевых значений частот собственных колебаний. Проведена коррекция свободной КЭ-модели ДПМ самолета Ту-204 по шести наборам частот собственных колебаний, определенных экспериментально. Целевые частоты при коррекции по каждому из шести наборов были достигнуты с высокой степенью точности.

About the authors

D. A Krasnorutskiy

Novosibirsk State Technical University; Siberian Aeronautical Research Institute named after S. A. Chaplygin

P. A Lakiza

Novosibirsk State Technical University; Siberian Aeronautical Research Institute named after S. A. Chaplygin

V. A Berns

Novosibirsk State Technical University; Siberian Aeronautical Research Institute named after S. A. Chaplygin

E. P Zhukov

Novosibirsk State Technical University; Siberian Aeronautical Research Institute named after S. A. Chaplygin

References

Statistics

Views

Abstract - 371

PDF (Russian) - 275