Применение аналога интеграла Коши-Лагранжа для расчета давления при моделировании трехмерного отрывного обтекания тел идеальной жидкостью Текст научной статьи по специальности «Физика»

УДК 532.5 + 517.9

ПРИМЕНЕНИЕ АНАЛОГА ИНТЕГРАЛА КОШИ-ЛАГРАНЖА ДЛЯ РАСЧЕТА ДАВЛЕНИЯ ПРИ МОДЕЛИРОВАНИИ ТРЕХМЕРНОГО ОТРЫВНОГО ОБТЕКАНИЯ ТЕЛ ИДЕАЛЬНОЙ ЖИДКОСТЬЮ

И. САТУФ, А.В. СЕТУХА

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

Ключевые слова: несжимаемая жидкость, вихревые методы, отрывные течения.

В работах С.М. Белоцерковского и его последователей [1 - 3] развиты численные методы моделирования трехмерного отрывного обтекания тел идеальной несжимаемой жидкостью, основанные на аппроксимации поверхностей тел и вихревых следов системами связанных вихревых отрезков и замкнутых вихревых рамок. При этом для вычисления давления используется интеграл Коши-Лагранжа.

Следует заметить, что при моделировании вихревого следа вихревыми рамками может происходить сильное растяжение сторон этих рамок. При этом вихревой след теряет возможность огибать препятствия, теряется возможность моделирования процессов самоорганизации вихревых структур. В работах [4, 5] описан подход, основанный на комбинированном использовании замкнутых вихревых рамок и изолированных вихревых элементов, который позволяет преодолеть указанные сложности. Однако использование интеграла Коши-Лагранжа для расчета давления в новом подходе оказалось некорректным из-за наличия завихренности в потоке жидкости. В работе Дынниковой Г.Я. [6] был предложен аналог интеграла Коши-Лагранжа и выведены аналитические формулы для расчета давления в нестационарном вихревом течении идеальной несжимаемой жидкости. В настоящей работе предложены формулы для расчета давлений на поверхности тел в рамках указанного выше комбинированного метода, полученные на основе результатов Дынниковой Г.Я.

1. Метод изолированных вихревых отрезков

Рассматривается нестационарное отрывное обтекание тела или системы тел в рамках модели идеальной несжимаемой жидкости.

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

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

ними сторонами таких рамок, через некоторое заданное число шагов после своего образования на линии отрыва преобразуется в изолированный отрезок, длина которого не должна увеличиваться. В результате отрезки, связанные в замкнутые рамки, разрываются, и образуется участок вихревого следа, моделируемый изолированными вихревыми отрезками. Численная схема такого метода описана в работах [4, 5].

Рис. 1

В данной статье рассмотрен вопрос о вычислении давления внутри жидкости и на поверхности обтекаемых тел при моделировании течения данным методом. При этом будем исходить из гипотезы, что дальняя зона вихревого следа, аппроксимируемая изолированными вихревыми отрезками, представляет собой область с объемным распределением завихренности, течение в которой подчинено нестационарным уравнениям движения идеальной несжимаемой жидкости.

2. Расчет нагрузок

Пусть Х1 есть поверхность обтекаемого тела (сумма поверхностей обтекаемых тел, если их несколько). Рассмотрим течение жидкости вне обтекаемого тела (системы тел), в котором имеется вихревая пелена (система вихревых пелен), на которой терпит разрыв касательная компонента вектора скорости, и одна или несколько ограниченных вихревых областей. Пусть Х2 -суммарная подвижная поверхность, на которой расположены вихревые пелены (ближний вихревой след); О - объединение всех вихревых областей (дальний вихревой след). Всюду в области течения вне системы областей О и вихревой пелены Х2 завихренность равна 0.

Будем считать, что в пространстве введена декартова система координат и каждую точку течения жидкости можно рассматривать как элемент х = (х1, х2, х3) пространства Я3. Предположим, что для поля скоростей жидкости w = '(х, 1) в каждый момент времени 1 справедливо следующее представление

' = + '1 + '2 + '3; (1) ' 1 =Уф1, ф1(х,1) = | - у)&(у,1)ёоу, 1 = 1,2, Б(х - у) = 1 1

Эп у

4р х - у

г 1 1

w3(x,t) = j w(y,t) x V(x - y)dy, V(x - y) = -VxF(x - y) = —-,,

W 4p lx - y|

w¥ - скорость жидкости на бесконечности; V = (Э/Эх1з Э/Эх2 ,Э/Эх3) - оператор Набла, интегралы по поверхностям понимаются как поверхностные первого рода, по области W - как тройной интеграл.

При этом векторные поля w1 и w2 являются градиентами потенциалов двойного слоя, размещенного на поверхностях обтекаемых тел S1 и вихревой пелене S2, как в модели обтекания тел с образованием вихревых пелен, описанной в работе [7], а w3 есть векторное поле, индуцируемое завихренностью w, сосредоточенной в дальнем вихревом следе. Будем считать, что функция ю определена во всем пространстве так, что w(x) = 0 вне системы областей W. При этом во всей области течения вне вихревой пелены S2 векторное поле w, определяемое формулой (1), удовлетворяет равенству rotw = ю.

Предположим также, что во всей области течения вне вихревой пелены S2 выполнены уравнения неразрывности и Эйлера (в форме Громеки-Ламба)

2

p w

divw = 0 , Эw/ Эt + w x w = -V(— + —), (2)

P 2

где p - давление; p - плотность жидкости.

Заметим, что при этом выполнено также уравнение переноса завихренности ЭюЭ = rot[wxw]. Из выражения (1) следует, что

Э^/ Э; = Э; + Эф 2/ Э;) + Эw 3/ Э;. (3)

Для производной Эw3/Эt в точке x в момент времени t справедливо выражение

Эw3_(x,t) = jЭЮ(У,;) x V(x - y)dy = j rot [w(y, t) x w(y,t)]x V(x - y)dy .

Э; W Э; W

Как показано в работе Дынниковой Г.Я. [6], для данного интеграла справедлива формула

j rot [w(y, t) x w(y, t)]x V(x - y)dy = Vl(x, t) + w(y, t) x w(y, t), (4)

W

где I(x, t) = j [w(y, t) x V(x - y)]w(y, t)dy .

W

Из выражений (2) - (4) получаем

V(p/ p + w 72 + Эф^ Э; + Эф2/ Э; +1) = 0. Так как в силу ограниченности области W выполнено условие I(x, t) ® 0 при |x| ® ¥, а потенциалы ф1 и ф2 также затухают на бесконечности, считая, что давление на бесконечности имеет предел p¥, из последнего равенства получаем следующий интеграл для нахождения давления

Р = + - -ЭФ1 -Ф -1 (5)

p p 2 2 Э; Э; . ( )

При численном моделировании течения описанным методом вихревых отрезков составляющие поля скоростей wi (i = 1, 2, 3) из правой части формулы (1) аппроксимируются в момент времени tk выражениями (рис. 1):

w1 = xgk(l)Wlme[l], w2 = xg(l)Wlme[l], (6)

leMj leM2(tk)

где M1 - множество вихревых рамок, моделирующих поверхность тела; M2(tk) - множество вихревых рамок, моделирующих поверхность ближней зоны вихревого следа в момент времени tk, сумма берется по всем рамкам l, моделирующим соответствующую поверхность; g(l) - циркуляция рамки l;

Wllne [l](x) = J dIy X V(X - y) l

у - элементарный вектор обхода кривой 1 - контура рамки; у - точки кривой 1;

' з = i Г(1ЖИпе [11

1еМз(1к)

где М3(1к) - множество вихревых отрезков, моделирующих дальнюю зону вихревого следа в момент времени 1к; Г(1) - циркуляция вихревого отрезка 1. При этом потенциалы ф1 и ф2 в момент времени 1к вычисляются по формулам:

jk =-Zgk(l)Uframe[о] jk =- Zg(l)Uframe[о], Uframe[о] = , (7)

leMj leM2(tk) о ЭП y

о - поверхность ячейки разбиения, охватываемая рамкой l.

Тогда производные по времени этих потенциалов можно вычислить по формулам:

j = -z ^^U_e [о]

Э" leM, At (8)

j = - egOf^-"J" zg(l)Uftame[0]

Э" leM2(tk-,) Э" DtleM2(tk)\M2(tk-,)

здесь * \ * - разность множеств (M2(tk)\M2(tk-1) - множество рамок, образовавшихся на последнем шаге интегрирования по времени); At - шаг интегрирования по времени (рис. 1).

3. Нахождение производной по времени от потенциала вихревой рамки

Пусть о = o(t) - подвижная ячейка, краем которой является замкнутая кривая l = l(t) . При этом контур l является подвижным. Пусть в каждый момент времени t этот контур задан параметрически x = y (s, t), где s - некоторый параметр. При этом w(s, t) = 9y(s, t)/9t есть скорость движения точки контура, соответствующей параметру s. Рассмотрим потенциал двойного слоя

j(x,t) = j af(p)dо.

о(") Эпу

При этом по закону Био-Савара справедлива формула

Vj(x, t) = - j t(s, t) X V(x - y(s, t))ds, (9)

l(t)

где x - точка вне поверхности о; T(s,t) = 9y(s,t)/9s - направляющий вектор на контуре l (здесь мы предполагаем параметр s введенным так, что вектор t(s, t) направлен в направлении положительного обхода контура l). В этом пункте мы докажем следующую формулу

j(x,t) = J[t(s,t)XV(x-y(s,t))]w(s,t)ds, xe R3\(aul). (10)

Э" l(t)

Для доказательства формулы (10) сначала продифференцируем выражение (9) по времени под знаком интеграла

V Мм) = - j Эт(М) X V(x - y(s, t))ds - j T(s,t) x 9V(x - y(s,t)) ds.

Э" l(t) Э" l(t) Э"

Заметим, что Эт/Э" = 3w/Эs, ЭV(x - y)/Э" =-(wVx )V(x - y), ЭV(x - y)/Эs =-(tvx )V(x - y), где оператор V x есть оператор Набла, вычисляемый по переменной x, y = y(s, t). Далее

j 3w/Э8 X Vds = j Э [w X V^sds - j w хЭ'У/Эsds, l(t) l(t) l(t)

где V = V (x - y(s,t)).

Так как первый из интегралов в правой части равен 0 в силу замкнутости контура 1(1), имеем

Эф

эТ

V ^ = J [т х (w Vx) V(x - y) - w x (тVx) V(x - y)]ds.

i(t)

Рассмотрим фиксированную точку y на контуре l в фиксированный момент времени. Пусть v(x) = V(x - y), J = т x (wV)v - w x (tV)v, w = (wj, w2, w3) , т = (tj, t2, t3). Тогда

3

J = z Witijuij, где u¡j = ej x (e¡V)v - e;V(ejV)v,

i,j=1

ei, i = 1,2,3 - орты координатных осей. Непосредственной проверкой устанавливаем, что выполнена формула uij =V([ej х v]ei), i,j = 1,2,3. Отсюда следует, что градиенты выражений в

правой и левой частях формулы (10) равны и, значит, эти выражения в каждый момент времени отличаются на константу. Для доказательства формулы (10) остается заметить, что оба указанных выражения стремятся к 0 на бесконечности и указанная константа равна 0.

4. Численное нахождение давления

Рассмотрим вопрос о приближенном вычислении давления в точке x в области течения на основании формулы (5). Скорость w вычисляется по формуле (1) с использованием выражений (6) и (7). Для нахождения производной Эф^ 9t предлагается использовать выражение (8). Если поверхность тела разбита на четырехугольные ячейки, то для вычисления потенциалов двойного слоя Uframe [о] каждая ячейка разбивается на два треугольника, вершины которых обходятся в том же порядке, что и в исходной рамке. Если поверхность S есть треугольник ABC, то потенциал Uframe [s] в точке x = M вычисляется как телесный угол, под которым видна поверхность S из точки M: Uframe[S](x) = q(a + b + g-p), где a,b,g - величины двугранных углов в

тетраэдре ABCM при ребрах AM, BM, CM соответственно, q = 1 при p > 0, q = 0 при

® ® ®

p = 0, q = -1 при p < 0 , где p = (ACx AB)MA .

Интеграл I приближенно вычисляется по формуле

I = z G(l)[r x V]w,

leM3(tk)

где r - вектор, соединяющий начало и конец вихревого отрезка; V = V(x - y), y - середина вихревого отрезка; w - скорость жидкости в точке y, вычисленная по формулам (1), (6), (7). В силу формулы (10) для производной Эф2/Эt справедливо приближенное выражение

ЭФ2М = z [r x V]w--1 z g(l)Uframe [O ,

leM2(tk-1) At leM2(tk)\leM2(tk-1)

где первая сумма берется по множеству M'2(tk-1) всех вихревых отрезков, из которых состоят вихревые рамки, существовавшие на поверхности s2 в момент времени tk-1 , причем, совпадающие отрезки соседних рамок при этом объединяются в единый отрезок, T(l) ; r и w - суммарные интенсивности, направляющие векторы, и скорости центров таких отрезков в момент времени tk (рис. 1). Для каждой вновь образовавшейся ячейки вихревого следа о, по контуру которой расположена вихревая рамка l е M2(tk)\M2(tk-1), можем записать

D Uframe №) = S = D V(x - y К = [w x r]V(x - y) = [r x V(x - y )]w ,

At At Эп y At

где Б - площадь данной ячейки; п0 - орт вектора нормали к этой ячейке; г - вектор, соединяющий начало и конец вихревого отрезка, противоположного стороне, лежащей на линии отрыва; у - середина этого отрезка; w - скорость жидкости в точке у . При этом мы использовали приближенное равенство Бп0 = w х г . Таким образом, суммируя выражения для интеграла I и производной Эф2/ Э1, получаем формулу

Эф2/ Э1 +1 = i Г(1)[г х У^, (11)

1еМ(1к)

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

Расчет краевых значений давления на поверхности обтекаемого тела производится по формуле (5) в контрольных точках (точках коллокации), тех же самых, которые используются для записи граничного условия непротекания поверхности тела. При этом в формуле (5) полагаем w = w0 ± А^2, где w0 - прямое значение вектора скорости жидкости, вычисленное по формулам (1), (6), (7); знак " +" берется, если давление считается на положительной стороне поверхности, " — " - в противоположном случае; Аw скачок скорости на поверхности обтекаемого тела Х1, который приближенно определяется по формуле [7]

Аw = у х п , у = i Г^/Б, где у - приближенная интенсивность вихревого слоя в данной точке; п - вектор нормали к рамке, на которой лежит рассматриваемая точка, сумма берется по всем сторонам вихревой рамки, в которой находится рассматриваемая точка; г - векторы, соединяющие начала и концы вихревых отрезков, лежащих на сторонах рамки, Г; = g — , где g - интенсивность данной рамки; - интенсивность соседней рамки, стыкующейся с рассматриваемой по вектору г;; Б - площадь рамки.

Производная Эф^ Э1 вычисляется по формуле

Эф1/ Э1 = (Эф1/ Э1)0 ± — gk—11)/2, где производная (Эф^ Э1;)0 вычислена по формуле (8), gk и gk—1 - интенсивность рамки, на которой лежит рассматриваемая точка в моменты времени 1;к и 1;к—1, знаки ± соответствуют положительной и отрицательной сторонам поверхности. Сумма Эф2/ Э1 +1 вычисляется по формуле (11) так же, как и для точки в потоке жидкости.

Наконец, заметим, что при вычислении разности давлений на тонкой поверхности слагаемое Эф2/ Э1 +1 сокращается в силу непрерывности, поэтому перепад давления можно вычислять по той же схеме, что и в методе вихревых рамок [7].

5. Примеры расчетов

На рис. 2 показаны результаты расчета аэродинамических характеристик прямоугольного крыла с удлинением 1 = 2,5 в диапазоне углов атаки от 0° до 90° при нулевом угле скольжения. Приведены зависимости от угла атаки коэффициентов нормальной силы Су и момента тангажа

ш2, которые вводились по формулам Су = 2Y/(pw^Б), ш2 = 2Мг/pw^Ь, где У - нормальная сила (направленная по нормали к поверхности крыла); М2 - момент тангажа, вычисляемый относительно оси, проходящей по передней кромке крыла; Б - площадь крыла; Ь - хорда крыла; р - плотность воздуха; w ¥ - как и ранее, скорость невозмущенного потока на бесконечности.

Для сравнения приведены данные эксперимента, полученные при продувках соответствующего крыла в аэродинамической трубе [8].

(точки - расчет; сплошная линия - эксперимент [7])

Рис. 2

Также был произведен расчет обтекания и распределений давления по каждой стороне квадратной пластины потоком, набегающим по нормали к пластине. На рис. 3 приведены распределения коэффициента давления р = р/pw¥ вдоль расчетной линии, проходящей по оси симметрии пластины, по стороне пластины, обращенной к набегающему потоку - р+ и по противоположной стороне - р . Также на рисунке приведены средние значения коэффициента давления на каждой стороне пластины, полученные в расчете, в сравнении с экспериментальными значениями [9].

1.0 ---■— —■---в

к ^^ / 0.5

/ \

-0.4 -0.3 -0.2 -0.1 0 / -0.5 0.1 0.2 0.3 0.4 ^

/ К / ■ ■—■-■— —■-" ■ \

!—■--■ "--вО

расчетная лнння

-0.5

0.5 я

средние значения коэффициента давления

расчет эксперимент [9]

Г 0.39 0.38

Т -0.82 -0.80

Рис. 3

Выводы

Приведенные результаты расчетов свидетельствуют об удовлетворительном совпадении с экспериментальными данными как суммарных нагрузок, так и распределений давления по поверхности тел, получаемых по предложенным формулам.

ЛИТЕРАТУРА

1. Белоцерковский С.М., Ништ М.И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. - М.:Наука, 1978.

2. Апаринов В.А., Дворак А.В. Метод дискретных вихрей с замкнутыми вихревыми рамками // Труды ВВИА им. проф. Н.Е.Жуковского. - 1986. - Вып. 1313. - С. 424 - 432.

3. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. - М.: ТОО "Янус", 1995.

4. Кирякин В.Ю. Моделирование обтекания объектов методом дискретных вихрей с представлением вихревой пелены изолированными вихревыми частицами // Научный Вестник МГТУ ГА, серия Аэромеханика и прочность. - 2008. - № 125 (1). - С. 78 - 82.

5. Апаринов А.А., Сетуха А.В. О применении метода мозаично-скелетонных аппроксимаций матриц при моделировании трехмерных вихревых течений вихревыми отрезками. // ЖВМ и МФ. - 2010. - Т. 50. - № 5. - С. 937 - 948.

6. Дынникова Г.Я. Аналог интегралов Бернулли и Коши-Лагранжа для нестационарного вихревого течения идеальной несжимаемой жидкости // Изв. РАН МЖГ. - 2000. - № 1. - С. 31 - 41.

7. Гутников В.А, Лифанов И.К., Сетуха А.В. О моделировании аэродинамики зданий и сооружений методом замкнутых вихревых рамок // Изв. РАН МЖГ. - 2006. - № 4. - С. 78 - 92.

8. Петров Е.Г., Табачников В.Г. Экспериментальное исследование аэродинамических характеристик прямоугольных пластин различного удлинения в широком диапазоне углов атаки // Труды ЦАГИ. 1974. - Вып. 1621. - C. 102 - 109.

9. Савицкий Г.А. Ветровая нагрузка на сооружения. - М.: 1972.

APPLYING THE CAUCY-LAGRANGE INTEGRAL ANALOGUE TO CALCUATING PRESSUSURE WHEN MODELING 3D SEPARATE FLOWS OF IDEAL FLUID PAST BODIES

Satuf I., Setukha A.V.

The approach presented which allows calculating pressure in fluid flow and at body surfaces when modeling 3D separate flows past bodies with the composed method of closed vortex frames and isolated vortex segments. Analytical formulas for calculating pressure, numerical method and computing results for sample problems are given.

Key words: incompressible fluid, vortex methods, separated flows.

Сведения об авторах

Сатуф Ибрагим, 1968 г.р., окончил КВВАИУ (1991), адъюнкт кафедры высшей математики ВУНЦ ВВС "ВВА им. проф. Н.Е. Жуковского и Ю.А.Гагарина", автор 3 научных работ, область научных интересов - вычислительная аэродинамика, вихревые методы.

Сетуха Алексей Викторович, 1966 г.р., окончил МГУ им. М.В.Ломоносова (1988), доктор физико-математических наук, профессор, заведующий кафедрой высшей математики ВУНЦ ВВС "ВВА им. проф. Н.Е. Жуковского и Ю.А.Гагарина", автор более 80 научных работ, область научных интересов - вычислительная аэродинамика, вихревые методы, интегральные уравнения, краевые задачи.