Диссертация (1025494), страница 9
Текст из файла (страница 9)
Для определения обобщенных координат {φ} оказываетсявозможным на каждом шаге интегрирования «припасовыванием» получатьзначения обобщенных координат с использованием аналитического решения ввиде [98]:φk (t i+1 ) =gkω2k+exp(−nωk ∆t)ω2k √1−n2{[ω2k √1 − n2 φk (t i ) − g k ] cos (ωk √1 − n2 ∆t)+ [ ωk φ̇k (t i ) + n(ω2k φk (t i ) − g k )] sin (ωk √1 − n2 ∆t)}φ̇k (t i+1 ) =exp(−nωk ∆t)n2{[√1 − n2 φk (t i ) − g k ] cos (ωk √1 − n2 ∆t)√1 −gk+ [ − (nφ̇k (t i ) + ωk φk (t i ))] sin (ωk √1 − n2 ∆t)}ωk(2.11)49Нелинейный позиционный характер нагрузки учитывается при этомвследствие изменения граничных условий (2.7) для уравнений гидродинамики, чтоприводит к изменению вектора правой части в (2.9).
Вследствие учета процессоввихреобразования, приводящих к непотенциальному характеру течения среды,изменение нагрузки нелинейно зависит от изменения положения, и скорости точекповерхности тела. Это делает невозможным разделение вектора {g} на отдельныеслагаемые, зависящие только от {φ̈}, {φ̇}, {φ}, с введением в рассмотрение матрицприсоединенных масс, аэродинамического демпфирования и аэродинамическихпроизводных как это обычно делается при решении задач аэроупругости. Вразработанной методике вектор {g} рассчитывается на каждом шаге сиспользованием метода вихревых элементов.Для формирования системы (2.3) и расчета перемещений узлов обтекаемоготела строится конечно-элементная модель.
Поверхность тела разбивается на конечных элементов.Одна из ключевых особенностей разработанной методики заключается вуниверсальности схемы формирования конечно-элементной модели. Сетка дляМКЭ может разрабатываться на доступных программных комплексах, которыесебя уже хорошо зарекомендовали в аэрокосмической отрасли, это MSC Nastran,ANSYS или с использованием программ с открытым исходным кодом, напримерCode Aster. Модель может быть достаточно сложной и содержать различные видыэлементов. Ограничением конечно-элементной модели является необходимостьзадания обтекаемой части конструкции в виде непрерывной замкнутойповерхности при помощи оболочечных или твердотельных элементов.
Балочные иодномерные элементы могут быть использованы в модели, однако их обтекание небудет рассчитано.2.3. Описание аэродинамической подсистемыИзвестны два подхода к изучению течения среды. В случае Эйлерова подходаконтролируются параметры среды в фиксированных точках пространства. Поток«протекает» через эти точки. В случае Лагранжева подхода контролируются50параметры среды в фиксированных точках потока.
Эти точки, жидкие частицы –маркеры перемещаются в НСК .В настоящей диссертации для решения задачи (2.1), (2.2), (2.5)-(2.7)применяется Лагранжев подход к изучению потока среды, поскольку он вместоиспользования сетки позволяет осуществить переход от системы уравнений вчастных производных к системе обыкновенных дифференциальных уравнений дляпараметров частиц, переносимых потоком среды.
Метод, в котором в качествечастиц применяются вихревые элементы - метод вихревых элементов (МВЭ)достаточно хорошо развит в настоящее время. Например, он подробно описан в[116].По теореме Гельмгольца о разложении векторных полей [82], полнаяскорость течения в любой точке представляется в виде⃗ =⃗ + ⃗ ,(2.12)⃗ = 2 – поле скоростей потенциального течения со скалярнымгде ⃗ = × ⃗⃗ – соленоидальное поле скоростей с векторнымпотенциалом , ⃗⃗ ().потенциалом Поскольку при решении задачи аэроупругости необходимо учитывать такиесложные гидродинамические эффекты, как вихреобразование и отрыв потока сгладкой поверхности, требуется рассмотрение сложных вихревых следов,геометрия которых не может быть сведена к тонким пеленам.
В этом случаенеобходимо представить поле скоростей в виде⃗ =⃗∞ + ⃗ ,⃗ определяется полем завихренностигде поле скоростей (2.13)⃗ (, ) =⃗ (, ).Применение операции ротора к обеим частям уравнения (2.2), послепреобразований и учета уравнения несжимаемости среды (2.1) [45] дает уравнениеНавье-Стокса, записанное в форме Гельмгольца [49]51⃗⃗ = (⃗ ∙ ∇)⃗,⃗ ∙ ∇)⃗ + ∇2 + ((2.14)Поскольку левая часть (2.14) есть субстанциональная производнаязавихренности, которая переносится жидкой частицей с маркером по полю⃗ , то (2.14) можно представить в лагранжевой формескоростей ⃗=⃗⃗ ∙ ∇)⃗,⃗ + ∇2 = ((2.15)Как видно из правой части (2.15), эволюция завихренности в области теченияскладывается из эволюции завихренности и ее диффузии, вызванной вязкостьюсреды.На обтекаемой поверхности с учетом граничного условия прилипанияуравнение (2.14) трансформируется в уравнение типа диффузии [47]⃗⃗ )| ,= ∇(∇[ ]∈ ∈(2.16)которое показывает, что среда, смежная со стенкой, вследствие прилипанияприобретает завихренность.
Эволюция завихренности на поверхности можетбыть представлена как процесс генерации новой завихренности, которыйхарактеризуется потоком завихренности с поверхности⃗ =Ω = ∇⃗,⃗(2.17)Таким образом, вязкость среды в уравнениях движения среды, во-первых,вносит вклад в эволюцию завихренности, и, во-вторых, определяет потокзавихренности при ее генерации вблизи поверхности тела. Принятое врассматриваемой задаче допущение о малой вязкости среды позволяет пренебречьвлиянием вязкости при эволюции завихренности в потоке - обнулить последнееслагаемое во втором уравнении (2.15).
Однако при этом поток завихренности (2.17)не обнуляется, что соответствует известному подходу Прандтля.52Наличие потока завихренности с обтекаемой поверхности обеспечиваеткинематическую согласованность полей скоростей и завихренности. Различныеподходы к вычислению потока завихренности подробно рассмотрены в шестойглаве монографии [47]. В настоящее время наибольшее распространение ввихревых методах получил подход Лайтхилла [60, 61]. Согласно этому подходу,генерация завихренности рассматривается как последовательность двух процессов:возникновения присоединенной к обтекаемой поверхности тонкой вихревойпелены и последующая диффузия завихренности из вихревой пелены в областьтечения, что обеспечивает для замкнутого тела удовлетворение условийприлипания [62].
Интенсивность вихревой пелены γ может быть найдена изусловия непротекания путем решения интегрального уравнения для точек наповерхности тела ∈ ⃗⃗⃗⃗ ( ) − ⃗ ( )],∫ ⃗ ∙ ( ( − ) × γ⃗ ) = ⃗ ∙ [(2.18)где ⃗ – нормаль к поверхности , ( − ) = [2( − 1)]−1 ( − )| − |−- градиент скалярной функции Грина для пространства размерности , - радиусвектор точки, по которой ведется интегрирование. Зависимость потоказавихренности от интенсивности вихревой пелены в (2.17) определяется какправило эмпирической моделью, пример которой можно найти в [63].Поскольку поле скорости (2.13) в области течения несжимаемой среды можетбыть восстановлено по известному полю завихренности с использованием законаБио-Савара⃗ ) ,⃗ (, ) = ⃗∞ + ∫ ( ( − ) × (2.19)описание течения несжимаемой среды в лагранжевой вихревой постановке (2.15)содержит только одну неизвестную функцию – поле завихренности, что упрощаетописание течения по сравнению с системой (2.1), (2.2).53Восстановление поля давления производится на основе известного полязавихренности.
Применение операции дивергенции к обеим частям уравнения(2.2), после преобразований дает уравнение Пуассона [65]⃗)⃗ ×∇2 = ∇ ∙ (V(2.20)для неизвестной величины(, ) = [p(, ) 1p122⃗ (, )| ] − [ ∞ + |⃗∞ | ]+ |ρ∞2ρ∞ 2(2.21)Уравнение (2.20) устанавливает связь между полным давлением в области⃗ [37].⃗ ×течения и так называемой массовой вихревой силой по Прандтлю Ω = VРешение (2.20) позволяет восстановить поле давления по известному полюзавихренности с использованием (2.17) для вычисления скорости. Формула длявычисления давления представляет собой аналог интеграла Коши-Лагранжа длязавихренных течений [9, 64].1122⃗∞ | − |⃗ (, )| + + Ω ](, ) = p∞ + ρ∞ [ |22(2.22)где – слагаемое, определяющее вклад от завихренности, рождающейся наповерхности, Ω - слагаемое, определяющее вклад от завихренности, находящейсяв области течения.
Для расчета коэффициентов , в случае плоскопараллельногои пространственного обтекания используются разные зависимости, которые будутуказаны в соответствующих разделах далее.Для перехода от системы уравнений в частных производных (2.15), (2.16) ксистеме обыкновенных дифференциальных уравнений используется методвихревых элементов (МВЭ).Поле завихренности представляется в виде суперпозиции элементарныхполей завихренности, называемых вихревыми элементами.⃗ (, ) = ∑ ⃗ ,=1(2.23)54⃗ =Вихревой элемент (ВЭ) представляет собой функцию трех параметров ⃗ , i ), где – маркер ВЭ – радиус-вектор характерной точки ВЭ в НСК; ℎ⃗ –( , ℎвектор ВЭ, связанный с маркером; – интенсивность ВЭ.Формулу (2.19) для расчета поля скорости по известным параметрам ВЭможно представить в виде суперпозиции скоростей, индуцируемых вихревымиэлементами⃗ ()),⃗ (, ) = ⃗∞ + ∑ () (, (), ℎ(2.24)=1⃗ ()) – скорость, индуцируемая ВЭ единичной интенсивностигде (, (), ℎпо закону Био-Савара. Конкретный вид этой функции определяется типомиспользуемого ВЭ и будет представлен в последующих главах.Принятые аппроксимации поля завихренности и скорости позволяютперейти от (2.15) к системе уравнений изменения параметров ВЭ [47]⃗ ( , )=⃗ℎ⃗ , ( = 1 … )= ℎ{ = 0(2.25)где – тензор деформации ВЭ.В силу сделанных преобразований в МВЭ граничное условие «набесконечности» (2.6) удовлетворяется тождественно, что делает данный методочень эффективным для расчета течений в безграничной области.Для удовлетворения граничного условия (2.7) на поверхности тела на каждомшаге интегрирования (2.25) проводится процедура генерации завихренности,состоящая из двух операций: решения уравнения (2.19) и формирования новых ВЭна всей поверхности тела .Для решения (2.19) используется панельный метод.
На поверхности строится сетка из Р плоских панелей для которых известно выражение единичного55⃗ ), ( = 1 … ). Условие непротекания, в заданных наскоса потока ∗ ( , , ℎпанелях контрольных точках записывается для неизвестных интенсивностейпанелей , ( = 1, … , Р ) в виде системы линейных алгебраических уравненийсовместно с уравнением баланса циркуляции.11⋮Р1[ 1………1Р⋮РР11 1⃗ 1⋮⋮⋮=,⃗ 1 РР0] { } { 0 }(2.26)⃗ ); – регуляризирующаягде = 1, … , Р , = 1, … , Р , = ∗ ( , , ℎ⃗ ))⃗ ,.⃗ = (⃗ − ⃗∞ − ∑V (, , ℎпеременная, =1С использованием найденных и геометрических параметров панелейопределяются параметры новых ВЭ:⃗ , = 1, … , b , , , ℎ(2.27)После формирования ВЭ они пополняют вихревую пелену - множество ВЭ вобласти течения. Таким образом на каждом шаге интегрирования размерностьсистемы (2.25) увеличивается: (t + ∆t) = V () + b .Начальными условиями для новых уравнений системы (2.26) являютсязначения (2.27). В нулевой момент времени начальным условием является полескоростей, реализующее бесциркуляционное обтекание тела за счет исходногораспределения интенсивностей по панелям обтекаемого тела.После удовлетворения граничных условий на шаге интегрирования могутбыть определены коэффициенты , и восстановлено поле давлений по формуле(2.22).Для определения вектора нестационарных аэродинамических сил {G(р)} поизвестному давлению среды ( , ) находятся аэродинамические силы ад ,действующие на узлы конечно-элементной модели.