Диссертация (1025494), страница 10
Текст из файла (страница 10)
При этом разбиениеповерхности обтекаемого тела на панели и конечные элементы должно быть56согласовано: радиус-векторы узлов панелей , = 1, … , , совпадают с радиусвекторами узлов конечно-элементной модели.2.4. Алгоритм решения задачи аэроупругостиРасчет переходного режима от начального до конечного момента временипроизводится в цикле из = / ∆ шагов.Алгоритм решения задачи, блок-схема которого представлена на Рис. 2.1, наодном шаге расчета при переходе от момента времени к моменту времени +1 = + ∆, состоит из следующих основных операций:1.Расчет параметров панелей обтекаемого тела в деформированномсостоянии на основе вектора перемещений {}.2.Построение расчетной схемы на поверхности обтекаемого тела иформирование СЛАУ (2.26)3.РешениеСЛАУ(2.26)формированиеначальныхпараметроврожденных ВЭ (2.27)4.Пополнение вихревой пелены ВЭ, рожденными на шаге расчета.Увеличение размерности системы (2.25). Новые ВЭ при генерации отодвигаютсяот панелей по нормали на малое заданное расстояние 5.Определение поля давления по формуле (2.22) и расчет нагрузок вузлах конечно-элементной модели6.Интегрирование системы дифференциальных уравнений динамикиупругой системы (2.3).7.Интегрирование системы дифференциальных уравнений движения ВЭ(2.25) методом Эйлера первого порядка с постоянным шагом ∆.В ходе этой операции полученные новые параметры ВЭ корректируются:- исключается попадание маркеров ВЭ внутрь обтекаемого тела за счетошибок округления.- не учитывается удлинение векторов ВЭ, которое превышает заданнуювеличину ∆ .57Рис.
2.1. Блок схема алгоритмаПосле интегрирования системы (2.25) проводится реструктуризациявихревой пелены, направленная на ускорение расчетов путем сокращения числа ВЭ[66]:- объединяются ВЭ, дистанция между маркерами которых меньше заданнойвеличины и косинус угла между векторами которых не меньше заданнойвеличины .-впеленесохраняютсятолькоВЭ,имеющиеминимальнуюрассматриваемую интенсивность. Элементы, интенсивность которых меньшезаданного порогового значения , удаляются- удаляются ВЭ на удалении от обтекаемого тела, превышающем заданнуювеличину ; данная процедура часто используется в МВЭ, поскольку в формуле(2.22) вклад в давление вихревого влияния элементов расположенных вдали от теластановится пренебрежимо мал по сравнению с близко расположенными ирожденными на шаге элементами.Программная реализация алгоритма на языке C++ в программе MDVDDпозволяет проводить параллельные вычисления по технологии MPI.582.5. Выводы по главе 2С учетом сделанных допущений, описанная в главе 2 методика решениясвязанной задачи аэрогидроупругости отличается тем, что позволяет:1.Рассчитывать малые колебания упругого обтекаемого тела методомразложения по собственным формам колебаний.
При этом формы колебанийдолжны быть получены на основе конечно-элементной модели любой сложности,в которой граница обтекаемого тела задана в виде непрерывной поверхности.Конечно-элементная модель может быть получена в любом известном пакете CAE.2.Рассчитыватьобтеканиедеформируемоготелаиопределятьаэродинамические нагрузки лагранжевым бессеточным методом вихревыхэлементов, что позволяет решать задачу обтекания в безграничной расчетнойобласти.3.Моделировать взаимодействие тела с потоком в случае, когдаплотность конструктивных элементов существенно больше плотности набегающейсреды (воздуха).
Взаимодействие определяется потоком завихренности со всейповерхности обтекаемого тела, который зависит от перемещения и скоростидвижения обтекаемой поверхности.4.Уменьшить время расчета за счет того, что вычислительные ресурсыиспользуются для расчета динамики вихревых элементов в спутном следе, числокоторых существенно меньше, чем число ячеек, необходимое для расчета задачисеточным методом.59Глава 3. Решение тестовых двумерных задач аэроупругостиВданнойглаверассматриваетсяплоскопараллельноеобтеканиеупругодеформируемых профилей: кругового и эллиптического. Производитсяверификацияразработанногопрограммногокомплексанаизвестныхэкспериментальных данных для абсолютно жестких профилей.
После этогопроизводится моделирование обтекания профилей, с учетом их податливостивнешней аэродинамической нагрузке.3.1. Расчетная схема упругой подсистемыМоделирование плоскопараллельного обтекания деформируемых профилейпроизводится при следующих допущениях:1.Упругодеформируемый профиль разбивается на участков.Распределение узлов является равномерным по всей длине контура, как показанона Рис. 3.1.321Рис. 3.1. Конечно-элементная модель профиля2.Каждый участок профиля представляет собой балочный конечныйэлемент, имеющий 3 степени свободы в каждом узле: растяжение , прогиб ,угол поворота и обладающий постоянной изгибной жесткостью по длине ,60равномерно распределенной постоянной погонной массой и длиной (Рис. 3.2).Рис.
3.2. Балочный конечный элементПоперечное сечение балочного конечного элемента является квадратом состороной (Рис. 3.3).Рис. 3.3. Сечение балочного конечного элемента3.Вводится неподвижная система координат (НСК) , начало которойсоответствует центру масс (ЦМ) профиля (Рис. 3.1). Движение набегающего потока⃗∞ = { , } и постоянной плотностью ∞ задается вс постоянной скоростью НСК.Вводится локальная система координат (ЛСК) ′′′ с началом в узлеконечного элемента. Координаты точки ′ (начала отсчета ЛСК) в НСК задаетсярадиус-вектором = {′ , ′} , ось ′′ направлена вдоль срединной линиибалочного конечного элемента, ось ′′ перпендикулярна срединной линиибалочного конечного элемента. Поворот конечного элемента в каждом из узлов отЛСК к НСК задается матрицей поворота [59][ ] = [ 0−000]1(3.1)613.2. Расчетная схема метода вихревых элементовВЭ в случае плоскопараллельного течения представляет собой бесконечнотонкую вихревую нить, перпендикулярную плоскости течения.В качестве вихревого элемента используется вихрь Рэнкина.
Это позволяетиспользоватьВЭ,которыйпредставляетсяввидевихревойнити,перпендикулярной плоскости течения, в виде круглых ВЭ фиксированного идостаточно малого радиуса.⃗ направленДля плоскопараллельного течения вектор завихренности перпендикулярно плоскости течения. В связи с этим, далее будем использовать⃗ (, ) на нормаль ⃗ к плоскости течения. Так как впроекцию (, ) вектора диссертационной работе рассматривается Лагранжев подход, то по теоремеТомсона завихренность в точке () не изменяется и () = (()) = .⃗ (, ) в области течения, используяПо известному полю завихренности закон Био-Савара [50], мы можем восстановить поле скорости в любой точке плоского течения:⃗ ( − ) × ⃗ ( ) d⃗ (, ) = ∫ (3.2)где⃗ ( − ) =1 ⃗ × ( − )2 | − |2(3.3)Так как при плоскопараллельном обтекании вектор ВЭ остается постоянным,⃗ = {0,0,1}ℎ⃗ℎ=0то система (2.25) вырождается в дифференциальное уравнение⃗ ( , )=(3.4)62Для численного решения уравнения (3.4) используется метод первогопорядка с временным шагом ∆.Суммарная завихренность в потоке среды в силу теоремы Томсона равняетсянулю.∑ = 0(3.5)=1Для численного моделирования процесса генерации завихренности строитсяследующая расчетная схема (Рис.
3.4).Рис. 3.4. Расчетная схема метода вихревых элементов для профиляВ качестве панелей профиля используются участки с узлами ( , +1 ). Накаждой панели выбирается контрольная точка на середине панели и точкарождения вихревого элемента с радиус-вектором ⃗⃗⃗⃗ = + ∙ ⃗ , где ⃗ - нормальк панели, - смещение от панели (Рис. 3.4). В контрольных точках контролируетсяусловие непротекания, а в точке рождения вихрей, для моделированиягенерируемой завихренности, устанавливаются ВЭ, интенсивности которыхподлежат определению.63Для вычисления давления в завихренных течениях жидкости применяетсяаналог интеграла Коши-Лагранжа [6]⃗∞2 ⃗2(, ) = ∞ + ∞ ( −+ − ),22(3.6)где⃗ ( − ) × ⃗ ( )) ∙ ⃗ ( , ) d = ∫ ((3.7) = ∫ (, )̇ ( ) , – слагаемое, учитывающее вклад в давление от движения завихренности вобласти течения; – слагаемое, учитывающее вклад в давление от рождениязавихренности.̇ ( ) – скорость рождения завихренности в точке ∈ , (, ) – потенциалскорости замкнутой вихревой линии единичной циркуляции, используемый прирасчете вклада в давление в точке с радиус-вектором от потока вихрей поповерхности панели.Сила, действующая на панель, направленная по нормали к ней, определяетсякак = − ∫ d(3.8)Для вычисления нестационарных аэродинамических нагрузок на i-ый узел⃗ , ) определяется в контрольных точкахконечно-элементной модели давление (⃗ ), = 1, … , П , и затем к каждому узлу прикладываетсяповерхности (аэродинамическая сила, распределенная по узлам панели.Аэродинамический момент относительно некоторой точки 0 равен⃗⃗ = − ∫( − 0 ) × (⃗ ) dгде ⃗ – единичная нормаль к панели профиля.(3.9)64Затем к каждому узлу прикладывается аэродинамический момент,распределенный по узлам панели.Сила лобового сопротивления и подъемная сила вычисляются какпроекции на оси НСК.
Аэродинамическим моментом является проекция вектора⃗⃗ , взятая с обратным знаком, на ось аппликат, направление которой соответствуеттому, что система координат правая.Длявычислениябезразмерныхаэродинамическихкоэффициентовприменяются следующие формулы:С () =2 ()2 ()2 ()(),С=,С()=−∞2∞2 2 ∞2(3.10)где – хорда профиля.Коэффициенты являются нестационарными. Для получения стационарныхзначений полученные зависимости будем осреднять по времени на достаточномпромежутке счета.3.3.
Результаты методических расчетовДля проведения верификации разработанного программного комплексарешения связанной задачи аэрогидроупругости был проведен ряд тестовых задач.Исследовалось обтекание профилей, для которых известны экспериментальныеданные стационарных аэродинамических коэффициентов лобового сопротивленияи подъемной силы.Для верификации метода конечного элемента были сопоставлены результатырасчета с результатами коммерческого пакета SolidWorks CosmosWorks.Моделированиеидальнейшееисследованиевсехтестовыхзадачпроводилось с использованием разработанной автором программного комплексаMKE2D.653.3.1.
Обтекание кругового профиляДля проверки адекватности разработанной методики и программногокомплекса были проведены серии расчетов, при которых профиль являлсяжестким, и далее, после сверки полученных результатов с экспериментальнымиданными были проведены расчеты, когда профиль являлся упругодеформируемым.Недеформируемый круговой профильВсе параметры расчетной схемы безразмерные и приведены в Таблице 2.Известно, что при обтекании кругового профиля в аэродинамическом следеза ним образуется так называемая дорожка Кармана. Для чисел Рейнольдса вдиапазоне = 103 … 105 безразмерная частота схода вихрей с цилиндра являетсяквазипостоянной и характеризуется числом Струхаля ℎ ≈ 0,2 , при этомустановившееся значение коэффициента лобового сопротивления ≈ 1,2, аосредненное установившееся значение коэффициента подъемной силы = 0.Число Струхаля вычисляется по формуле:ℎ =∙,∞(3.11)где – частота схода вихрей, – характерный размер тела, ∞ – скоростьсреды.Вид вихревой дорожки, образующейся за круговым профилем, в моментвремени показан на Рис.
3.5.66Рис. 3.5. Вихревой след за круговым профилемПри указанных параметрах расчетной схемы для жесткого цилиндра былиполучены установившиеся значения аэродинамических коэффициентов =1,196 и = 0. Амплитуда колебаний подъемной силы в расчетах составила = 0,9. Спектр пульсации гидродинамических сил – моногармонический сбезразмерными частотами ℎ = 0,21 и ℎ = 0,4. Эти данные хорошосогласуются с известными экспериментальными данными, упомянутыми ранее[57].На Рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
