Диссертация (1025494), страница 12
Текст из файла (страница 12)
Например, если параметры упругого профиля настроены в резонанс счастотой схода вихрей Кармана для абсолютно жесткого кругового профиля, то засчет изменения параметров вихреобразования возможна отстройка частотывынужденных колебаний от резонанса. В то же время у упругого профиля,изначально отстроенного от резонанса с частотой схода вихрей Кармана, можетвозникнуть режим автоколебаний.4.Расчеты показывают, что численное интегрирование уравненийдинамики упругой подсистемы требуют использовать малый шаг по времени длявысокочастотных вибраций.84Глава 4. Решение тестовых трехмерных задач аэроупругостиВданнойглаверассматриваетсяпространственноеобтеканиедеформируемых тел. Производится верификация разработанного программногокомплекса на известных экспериментальных данных для цилиндрическойоболочки.
После верификации производится решение ряда модельных задачобтекания конструкций, с учетом податливости поверхности тела внешнейаэродинамической нагрузке.4.1. Расчетная схема упругой подсистемыПроцесс построения конечно-элементной модели и получения глобальныхматриц масс [], демпфирования [] и жесткости [] сводится к следующимшагам:1.В пакете MSC PATRAN создается геометрическая модель обтекаемого2.Данные по модели выгружаются в файл с расширением .bdf3.Считываются координаты узлов сетки, что является координатамителавершин панелей МВЭ.4.Из пакета MSC NASTRAN определяются собственные формыколебаний оболочки с помощью решателя SOL 103. Поскольку матрица массы вMSC NASTRAN редуцируется к единичной, то в результате матрица жесткостипредставляет собой диагональную матрицу с квадратами собственных частот[2 ] .5.В результате, в файле с расширением .f06 содержатся матрицакоэффициентов собственных форм [], глобальные матрицы масс [],демпфирования [] и жесткости [], а также координаты узлов сетки оболочечныхконечных элементов.854.2.
Расчетная схема метода вихревых элементовВ отличие от плоскопараллельного течения, при моделировании эволюции⃗ () изменяются взавихренности в пространственном случае векторы ВЭ ℎсоответствии с уравнением движения среды (2.2).В качестве ВЭ выступает симметричный вортон-отрезок, который подробноописан в работе [101]. Структура поля завихренности ВЭ может быть получена каксуперпозиция полей точечных вортонов, составленных в прямолинейный отрезоктак, что векторы точечных вортонов направлены от начала к концу этого отрезка.После интегрирования поле завихренности симметричного вортона-отрезкапредставляется как суперпозиция полей:⃗)=⃗)+⃗)⃗ ∗ (, 0 , ℎ⃗ ∗0 (, 0 , ℎ⃗ ∗/ (, 0 , ℎ(4.1)⃗ ) - завихренность на отрезке⃗ ∗0 (, 0 , ℎгде ⃗ ), −ℎ < < ℎ,( − 0 − ℎ⃗ ) = ( − − ℎ⃗ ∗0 (, 0 , ℎ⃗ ), = −ℎ, = ℎ,02⃗ , 0 + ℎ⃗]{0, ∉ [0 − ℎ(4.2)⃗ ) - дополнительная завихренность⃗ ∗/ (, 0 , ℎ⃗)=×(⃗ ∗/ (, 0 , ℎ1112 ⃗⃗⃗⃗ ) =−()44 |1 |3 |2 |3(4.3)Поле (4.3) подобно векторному полю диполя, у которого источник и сток⃗разнесены на расстояние 2ℎ.
Семейство вихревых линий для ВЭ, у которого = и 0 = 0, ℎ = 1, показано на Рис. 4.1. Как видно из рисунка, поле завихренностиотрезка (4.2) замыкается линиями поля (4.3), что обеспечивает для данного ВЭвыполнение теорем Гельмгольца.86Рис. 4.1. Вихревые линии симметричного вортона-отрезкаДля моделирования потока завихренности со всей поверхности обтекаемоготела в предлагаемой расчетной схеме вихревые отрезки, составляющие рамку,заменяются вортонами. При этом, образуется вортонная рамка (Рис. 4.2), у которойчисло рождающихся вблизи панели ВЭ совпадает с числом сторон панели .Интенсивность всех ВЭ, рожденных вблизи одной панели, считается одинаковой иравной циркуляции рамки .Рис.
4.2. Вортонная рамка87Для удовлетворения граничных условий поверхность обтекаемого телааппроксимируется многогранником, состоящим из П плоских панелей. Каждая jая -угольная панель задается своими параметрами: радиус-векторами вершин внеподвижной системе координат , = (1, … , ); радиус-вектором контрольной⃗ ; внешней единичной нормалью к панели ⃗ ; площадью панели ,точки направляющим вектором ребра панели , длиной стороны панели .Вычисляются следующие параметры панелей:⃗ =1∑ ⃗ ==1 = − | − |1 × 2|1 × 2 | = | − | = (1, … , )(4.4) + 1,={1, < = Функция влияния ВЭ i- ой панели на j-ую панель определяется как⃗ , ) = ∑ ( , 0 , ℎ(4.5)=1Суммарное количество ВЭ, рождающихся на поверхности тела за один шаграсчета, равноП = ∑ (4.6)=1Для определения скорости, индуцируемой вортоном-отрезком используетсявыражение⃗ 0 (, ) =⃗ × 0 )(ℎ21⃗]−[()∙ℎ2||||21⃗ × 0 |4|ℎ(4.7)которое можно переписать в виде⃗ 0 (, ) =⃗ × 0 ; с = | |−2 [( 2 − 1 ) ∙ ℎ⃗]с ; = ℎ|2 | |1 |4(4.8)88Для выбора функции сглаживания особенности в (4.7) было исследованополе скоростей при приближении к оси вортона, заданной вектором .Исследование показало, что сглаживающая функция для вортона-отрезка должнабыть не сферически симметричной, а цилиндрически симметричной, т.е.
зависетьот расстояния до оси вортона и, возможно, координаты вдоль вортона. Длясглаживания особенности был введен линейный закон убывания скорости. Такойподход аналогичен сглаживанию скорости для вихря Рэнкина:1с , ≥ ,4⃗(, 0 , ℎ, , ) = {1 ∗ ∗с , < ,4 (4.9)⃗ × ( ∗ − 0 ),∗ = ℎс∗⃗⃗ ∗ − 0 + ℎ ∗ − 0 − ℎ∗ |−2⃗ ]|= −[()∙ℎ∗∗⃗⃗| − 0 + ℎ| | − 0 − ℎ |⃗0 ℎ⃗ = + ( − 1) (ℎ− 0 ) , = | × 0 |⃗ℎ2∗В уравнении эволюции вектора симметричного вортона-отрезка⃗ℎ⃗ ∙ )⃗ = [(0 ) ∙ ℎ⃗]⃗ = (⃗ (0 )) ∙ ℎ= (ℎградиент поля скоростей можно представить как тензор деформации ВЭ(4.10)⃗ (0 ) = [(0 )]в виде суммы:[(0 )] = ∑[(0 , 0 )](4.11)=1⃗ - изменение вектора i-го ВЭ под действием j-го ВЭТогда [(0 , 0 )] ∙ ℎ89[(0 , 0 )] =(Частные(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0(0 (0 )) ⁄0производныекомпонентоввектораскорости(4.12))⃗ 0 (0 )покомпонентам радиус-вектора маркера ВЭ 0 в (4.12) могут быть найденынепосредственно из выражения (4.8).При расчете градиента скорости необходимо вводить сглаживание наотрезке ВЭ в правой части уравнений (4.10).
Эта процедура производитсяаналогично (4.9) с использованием вычисленных в точке ∗ функций с∗ , ∗ ,.[(0 , 0 )], ≥ ,[(0 , 0 )] = {[(0 , 0 , ∗ , с∗ , ∗ )], < ,(4.13)Каждому i-му узлу расчетной панели модели ставятся в соответствие Ппанелей на поверхности тела:П⃗ , ) ⃗ ,ад = − ∑ ( = 1, … , (4.14)=1Для вычисления давления среды на j-ую панель поверхности обтекаемоготела используется аналог интеграла Коши-Лагранжа:⃗2 ⃗2⃗ , ) = ∞ + ∞ ( ∞ −(+ − ),22гдеПП⃗ , )) ∙ ⃗ , ) = ∑ ⃗ (⃗ = ∑ ( ∙ (=1=1П1 = − )) ∑ (⃗⃗⃗ ∙ (⃗⃗⃗⃗⃗⃗ × ⃗⃗⃗⃗⃗4=1(4.15)90В формуле (4.15) телесный угол , под которым видна панель из точки срадиус-вектором ⃗ определяется как [104]: == 2 ∑ (=1 ∙ (⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) × ⃗⃗⃗⃗⃗|⃗⃗⃗ ||⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ∙ ⃗⃗⃗⃗⃗⃗ )|⃗⃗⃗⃗⃗⃗⃗⃗ ∙ ⃗⃗⃗⃗⃗ )|⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ )|⃗⃗⃗ | || | + ( | + ( | + ( ∙ ⃗⃗⃗⃗⃗) + 1; < ⃗ − ⃗; ⃗⃗⃗⃗⃗⃗ = ⃗⃗⃗ = ⃗⃗⃗⃗ − ⃗; ⃗⃗⃗⃗⃗ = ⃗⃗⃗⃗⃗ − ⃗; = (1, … , ), = {1; = 4.3.
Результаты методических расчетовАлгоритм решения, в случае пространственного обтекания, отличается отплоскогопараллельного только решением упругой системы. Для определенияобобщенных координат {} используются аналитические выражения (2.11).Для проведения верификации разработанного программного комплексарешения связанной задачи аэроупругости был проведен ряд тестовых задач.Исследовалось обтекание тел, для которых известны экспериментальные данныеаэродинамических коэффициентов лобового сопротивления и подъемной силы.4.3.1.
Обтекание жесткой цилиндрической оболочкиРассматривается модельная задача для проверки адекватности разработаннойметодики и программного комплекса для трехмерных задач с целью дальнейшегоперехода к расчету обтекания прототипа легких РН.Вектор скорости невозмущенного потока сонаправлен с осью неподвижной системы координат и имеет модуль ∞ . Оболочка по торцам заделана,причем торцевые поверхности считаются абсолютно жесткими . Эти поверхностиобтекаются потоком.
Ось цилиндрической оболочки совпадает с осью OZнеподвижной системы координат. Параметры расчетной схемы даны в Таблице 9.91Таблица 9.Параметры расчетной схемыНаименование, размерностьОбозначениеЗначениеСкорость невозмущенного потока среды, м/с⃗∞Плотность среды, кг/м3∞1000∞105∆0,02Радиус сглаживания поля скорости ВЭ, с0,04Величина дистанции объединения ВЭ, м0,7Пороговое удлинение вортона, б/р∆1,50,9Смещение точки рождения ВЭ от панели, м0,02Дальность расчета вихревого следа, м410−6Статическое давление в невозмущенном потокесреды, ПаШаг интегрирования по времени уравненийаэродинамической подсистемы, сДопустимый косинус угла между векторамивортонов для объединения ВЭ, б/рМинимальная рассматриваемая интенсивностьВЭ, 1/сСхема проведения эксперимента приведена на Рис.
4.3. В ходе экспериментастальная цилиндрическая оболочка, размещенная на движущейся тележке,двигалась в воде [15].92Рис. 4.3. Цилиндрическая оболочка, установленная на тележке [15]Параметры и механические свойства материала оболочки в расчетной схемесогласно данным эксперимента равны: декремент колебаний = 0,05, модульупругости = 1,98 ∙ 1011 Па, коэффициент Пуассона = 0,31, плотностьматериала оболочки = 7920 кг/м3 , диаметр оболочки = 0,1156 м, длинаоболочки = 0,512 м, толщина оболочки = 0,0003 м.В коммерческом пакете MSC Patran была создана расчетная схема оболочкис общим количеством панелей равным 1400 (40 панелей в окружном направлениии 32 панели по образующей соответственно, а также 120 панелей на торцахоболочки).Расчеты проводились для двух значений модуля скорости набегающегопотокасреды:∞ = 0,9 м⁄с,чтосоответствоваловэксперименте[15]докритическому режиму обтекания оболочки и ∞ = 4,0 м⁄с, что соответствовало93в эксперименте [15] закритическому режиму обтекания оболочки.
Шаг временибыл равен в первом случае ∆ = 0,008 с, а во втором ∆ = 0,001 с.На Рис. 4.4 представлены полученные в расчете и осредненные по 800 шагаминтегрирования графики распределения коэффициента давления С по окружностицилиндра для рассмотренных скоростей набегающего потока.