Нестационарное обтекание профиля с угловой кромкой Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 5, № 3, 2000

НЕСТАЦИОНАРНОЕ ОБТЕКАНИЕ ПРОФИЛЯ С УГЛОВОЙ КРОМКОЙ*

Д. Н. Горелов, Ю.С. Смолин Институт информационных технологий и прикладной математики СО РАН, Омск, Россия e-mail: gorelov@iitam.omsk.net.ru, smolin@iitam.omsk.net.ru

An algorithm of solution for the nonlinear initial boundary value problem of unsteady incompressible flow around a profile with the angular trailing edge is suggested. This algorithm allows one to obtain the hydrodynamic characteristics near the trailing edge with high accuracy. It is shown that the exact realization of the boundary conditions on the trailing edge gives a discontinuity solution which in the limiting case of an infinitesimal angle does not coincide with the solution for a sharp trailing edge.

Нелинейная начально-краевая задача нестационарного обтекания профиля является одной из наиболее трудных в теории крыла. В общем случае циркуляция скорости вокруг профиля меняется с течением времени, что порождает вихревой след за профилем. Граничные условия на следе нелинейны, а его форма заранее неизвестна и меняется с течением времени. На некотором расстоянии от профиля вихревой след теряет устойчивость. Однако вблизи профиля как согласно теоретическим расчетам, так и в эксперименте наблюдаются устойчивые вихревые структуры. Это позволяет решать некорректную начальную задачу (в идеальной несжимаемой жидкости — задача Коши для уравнения Лапласа) в некотором конечном интервале вблизи начального момента времени.

Теоретическое исследование таких задач удается провести только для начальной стадии схода вихревого следа с профиля [2, 6, 12]. Что касается алгоритмов численного решения нелинейных начально-краевых задач, то их разработке посвящено большое число работ (см., например, обзор Сарпкайя [8]). Общим для всех алгоритмов является линеаризация задачи в малой окрестности рассматриваемого момента времени и моделирование вихревого следа системой дискретных вихрей. Соответствующие краевые задачи (для каждого момента времени) решаются обычно методом дискретных вихрей. Реже применяется метод панелей, хотя он позволяет решать задачу с большей точностью, особенно в окрестности задней кромки профиля.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-00093.

© Д.Н. Горелов, Ю.С. Смолин, 2000.

течения (см., например, [1, 3, 7, 12]). Вызвано это прежде всего тем, что процедура линеаризации, применяемая последовательно для ряда моментов времени, оправдана только для профилей с острой кромкой (в виде точки возврата). Нестационарное же течение вблизи угловой кромки имеет существенную нелинейность, которая связана как с изменением направления схода вихревого следа с профиля, так и с нелинейной зависимостью между циркуляцией скорости и интенсивностью сходящих вихрей.

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

1. Постановка задачи

Рассмотрим в декартовой системе координат 0xy нелинейную начально-краевую задачу о произвольном движении профиля с угловой задней кромкой в безграничном потоке идеальной несжимаемой жидкости. Пусть с момента времени t = 0 контур профиля L(t) начал двигаться со скоростью, определяемой вектором U(x,y,t), (x,y) £ L(t). В общем случае его движения циркуляция скорости r(t) вокруг профиля меняется с течением времени, что порождает вихревой след — линию контактного разрыва Lw (t) (рис. 1). Движение жидкости вне профиля и вихревого следа предполагаем потенциальным.

Задачу будем решать в комплексной плоскости z = x + iy, выбирая в качестве искомой функции комплексную скорость v(z,t) = vx(x,y,t) — ivy(x,y,t), которая при сделанных предположениях является аналитической функцией в области течения жидкости.

Функция v(z,t) должна удовлетворять граничным условиям непротекания жидкости через контур L(t), затухания возмущенных скоростей в бесконечно удаленной точке, непрерывности давления и нормальной составляющей скорости жидкости при переходе через вихревой след Lw (t). Для определения циркуляции скорости r(t) необходимо иметь дополнительное условие — постулат Кутта — Жуковского об ограниченности скорости жидкости в задней кромке. Это условие выполняется, если вихревой след сходит с задней кромки.

y U

Lw (t)

Рис. 1. Постановка задачи.

Начальные условия задаются тем, что при £ < 0 профиль неподвижен ("О = 0), течение стационарно и вихревой след ьт отсутствует.

Комплексную скорость будем искать в виде

1 Г 7, 1 [ (а,г)аа ( .

где — скорость жидкости в бесконечно удаленной точке, £(з,£) £ ¿(¿), (а,¿) € (¿), 7(в, ¿), 7ад(а, ¿) — интенсивности вихревых слоев, моделирующих профиль и вихревой след, 5, а — дуговые координаты точек £(в,£), (а, ¿).

Выражение (1.1) позволяет свести задачу к определению функций 7(в,£), 7^(а, ¿) и линии тангенциального разрыва ьт (¿). При этом выполняется граничное условие в бесконечно удаленной точке. Для выполнения условий на ьт (¿) достаточно потребовать, чтобы вихревой слой на ьт (¿) был свободным и перемещался вместе с жидкостью. Это имеет место, если комплексная координата (а, ¿) частицы жидкости, сошедшей с задней кромки профиля (точки гв) в некоторый момент времени ¿0 (0 < ¿о < ¿), определяется решением нелинейного дифференциального уравнения

^ = Фт ,£) (1.2)

с начальным условием (0,£о) = ¿в(¿о)- Здесь ¿0 является лагранжевой координатой рассматриваемой частицы жидкости, а , ¿) — скорость движения вихрей в следе, равная полусумме предельных значений комплексной скорости при подходе к точке

Интенсивность сходящих вихрей 7в (¿) для любого момента времени связана с циркуляцией скорости вокруг профиля Г(£) равенством

аГ

— = — 7в №в, 7в = 7т (0,£), (1.3)

аг

где №в (¿) — скорость схода свободных вихрей с задней кромки.

Сход вихрей может происходить только по касательной к одной из сторон профиля. При этом для профиля с угловой задней кромкой ($ = 0) другая сторона профиля не может гладко стыковаться с вихревым следом, образуя угловую точку, в которой относительная скорость жидкости должна быть равна нулю. В результате скорость схода вихрей становится равной wв = |7в |/2 и соотношение (1.3) принимает вид

- = —7В ^. (1.4)

аг 7в 2 1 ;

Изложенная постановка задачи хорошо известна. При этом в работах, посвященных решению задач нестационарного обтекания профиля с угловой кромкой, применяется соотношение (1.4), а не условие равенства нулю относительной скорости жидкости в угловой точке соответствующей линии тока, на основании которого оно получено [6-8, 12].

Отметим, что нелинейная зависимость (1.4) интенсивности сходящих вихрей 7в от циркуляции скорости вокруг профиля Г(£) делает принципиально невозможным линеаризацию задачи нестационарного обтекания профиля с угловой кромкой даже при малых возмущениях потока (для профиля с острой кромкой скорость схода вихрей в этом случае

мало отличается от скорости набегающего потока, что позволяет линеаризовать исходное соотношение (1.3)).

Рассматривая теперь выражение (1.1) в точках контура L(t) и применяя к интегралу типа Коши формулы Племеля, получим

U(z,t) = U(z,t) - 2Y(s,t)e-*(s,i), z g L(t), (1.5)

где интеграл по L(t) в (1.1) понимается в смысле главного значения по Коши, U = Ux — iUy, 0(s,t) — угол между касательной к L(t) и осью x, s — дуговая координата точки z g L(t) (см. рис. 1).

Умножая (1.5) на exp(i0(s,t)) и разделяя в полученном выражении вещественную и мнимую части, придем к двум интегральным уравнениям относительно функции y(s,t):

Im{e*(s,i)[v(z,t) — U(z,t)]} = 0, z g L(t), (1.6)

2Y(s,t) + Re {e*(s,i) [u(z,t) — U(z,t)]} = 0, z g L(t). (1.7)

Уравнение (1.6) определяет собой условие непротекания жидкости через контур L(t) (равенство нулю нормальной составляющей относительной скорости) и является сингулярным интегральным уравнением 1-го рода. Уравнение (1.7) налагает ограничения на касательную составляющую относительной скорости жидкости вдоль контура L(t) и является интегральным уравнением 2-го рода с ограниченным ядром.

Уравнения (1.6), (1.7) имеют параметрическую особенность, связанную с толщиной профиля. Эта особенность проявляется в том, что по мере утоньшения профиля уменьшается расстояние между соседними точками на верхней Li и нижней L2 сторонах профиля. В результате интегральные уравнения, записанные отдельно для верхней и нижней сторон контура L, в предельном случае профиля нулевой толщины становятся одинаковыми.

Параметрическая особенность создает серьезные трудности при численном решении интегральных уравнений (1.6), (1.7) для тонкого телесного профиля, а также в окрестности задней кромки для профиля произвольной толщины. В работе [5] предложено преобразование независимых интегральных уравнений 1-го и 2-го рода в систему двух совместных интегральных уравнений, не содержащих параметрической особенности. Применительно к уравнениям (1.6), (1.7) эта система имеет вид

Im {ei0(si'i)[v(zi ,t) — U(zi ,t)] — e*(s2,i) [U(z2, t) — U(z2,t)]} = 0, (1.8)

2[Yi(si,t) — Y2(s2,t)] + Re {ei0(si'i)[v(zi,i) — U(zbt)] — е*(«>4) [U(z2, t) — U(z2,t)]} = 0. (1.9)

Здесь z1, z2 — соседние точки на контурах L1, L2 соответственно; s1, s2 — дуговые координаты точек z1,z2; Yi, Y2 — искомые функции (yi = Y(si,t), Y2 = Y(s2,t)),

v(zr ,t) = U. + / / / ,r = 1, 2. (1.10)

2пг J zr — Z1(s,t) 2пг J zr — Z2(s,t) 2пг J z — zw(a, t)

Li L2 Lw

2. Метод решения

Для решения нелинейной начально-краевой задачи (1.1) - (1.5), (1.8) - (1.10) применим известную процедуру дискретизации по времени, позволяющую для ряда последовательных

значений времени ¿1, ..., ¿п получать систему интегральных уравнений (1.8), (1.9) с известными контурами ¿(¿т), (¿т.), т = 1, ..., п. Эти уравнения будем решать методом панелей с линейным распределением интенсивности вихревого слоя по каждому элементу контура Ь. Что касается вихревого следа , то его п-й элемент, примыкающий к задней кромке и образованный вихрями, сошедшими с профиля за время ДГ^ — ¿п ¿п—1, также моделируем панелью с линейным распределением вихрей. Другие же элементы следа, образованные вихрями, сошедшими с профиля за промежутки времени Д¿m = ¿т — ¿т-1, т = 1, ..., п — 1, заменим системой свободных дискретных вихрей Гтт, интенсивность которых равна суммарной интенсивности вихревого слоя на рассматриваемых элементах (с течением времени эти элементы деформируются, но суммарная интенсивность вихрей на них остается постоянной).

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

Для принятого моделирования контура и вихревого следа на каждом шаге по времени в качестве искомых выступают величины: т^ — значения 7(^¿) на стыках панелей (] = 0, ..., N, г = 1, 2; г = 1 соответствует верхней стороне профиля, а г = 2 — нижней), и +1 — значение интенсивности вихревого следа в конце его первого элемента (в точке +1).

Комплексная координата гм+1 определялась выражением

гм+1 = гв (¿п) + Дт , = ^в(¿«-1^«

(2.1)

где Дт — длина п-го элемента (панели) вихревого следа, — угол схода вихревого следа с профиля, отсчитываемый от оси х.

В результате для каждого момента времени ¿п интегральные уравнения (1.8), (1.9) сводятся к системе алгебраических уравнений вида

2 N

ЕЕтГ

г=1 ^=0

Ьу I г,

(1) 0к

Ь(г)

— Ьг I г,

(2) 0к

+1 (¿О*) — +1 (г

к = 1,

N.

—а (г0кМ +а (г,

(2)

(2.2)

(1 — ^) (тк-1 — + ^ (7Г — 7к

(2)

„(1)

(2)

2 N

+ ЕЕ ^

г=1 ^=0

а;-

0к

— а(Г) ( г(2)

+

+7м+1

ам+1 г,

(1)

ам+11 г

(2)

= —сг,

Д1)

+ с г,

(2)

к = 1,

N.

(2.3)

Здесь величины а(г) , Ь(г) и ам+1 ^0^), ЬМ+1 определяют касательную и

нормальную составляющие скоростей (в контрольных точках , р = 1, 2), индуцирован-

(г) ( (р)Н

ных линейным распределением вихрей, соответствующим значениям 7^ и 7м+1; с I ^ ), а I — касательная и нормальная составляющие известной части скорости жидкости в

точках г0рк), включая скорости, индуцированные системой свободных вихрей Гт1, ^ — коэффициенты, определяющие положение контрольных точек на панелях.

Г

тп-1;

1

1

2

Тестовые расчеты показали, что положение контрольных точек в системах уравнений (2.2) и (2.3) целесообразно выбирать разным. Поэтому в системе (2.2), соответствующей интегральному уравнению 1-го рода, контрольные точки выбирались в середине панелей (^ = 0, 5), а в системе (2.3), соответствующей уравнению 2-рода, — на 3/4 длины каждой панели (^ = 0, 75).

2Ж уравнений (2.2), (2.3) содержат 2М+3 неизвестные величины 7^ и 7м+1. В качестве дополнительных трех уравнений выбирались: условие непрерывности вихревого слоя в передней кромке профиля

^ = (2.4)

соотношение (1.3)

г(гп-1) - г(гп) = 1б Д^, = ^мг - тМ )/2, тв

7!!' + Т^

(2.5)

и уравнение

Г(^„-1) - Г(£„) = (тб + 7м+1) Дт/2,

(2.6)

которое является следствием теоремы Томсона о сохранении циркуляции скорости по замкнутому жидкому контуру и определения суммарной интенсивности вихрей на первом элементе вихревого следа.

В случае угловой кромки ($ = 0) необходимо требовать выполнение дополнительного условия равенства нулю относительной скорости жидкости в угловой точке на линии тока. С учетом перехода угловой точки с верхней стороны профиля на нижнюю в зависимости от знака ¿Г/^; это условие принимает вид

7'!1'

0 при ¿Г/^ < 0,

(2) тМ

0 при ¿Г/^ > 0.

(2.7)

Нелинейность условия (2.7) и разрывный характер поведения локальных характеристик в угловой кромке при смене знака ¿Г/^;, что было выяснено в результате численного эксперимента, потребовали внесения изменения в описанный выше алгоритм расчета. Условие (2.7) удалось выполнить путем замены линейного сплайна для интенсивности вихревого слоя 7(а) на N — панелях, примыкающих к задней кромке (точке ^в), на квадратичный сплайн специального вида. А именно,

7 (а)

(г)

тМ -1

1 -

а

Д(г) Д

, (г)

+ тМ

к

а

1г

Д(г) Д

+ к

а

2г

Д(г) Д

1-

а

Д(г) Д

а е

0, Д

1, 2.

(2.8)

Здесь коэффициенты к11 = 1, к12 = 0, к21 = 0, к22 = 1 при > 0 и к11 = 0, к12 = 1, к21 = 1, к22 = 0 при < 0. В случае острой кромки ($ = 0) коэффициенты к1г = 1,

к2г = 0.

Отметим, что уравнения (2.5), (2.7) приводят к конечно-разностному аналогу соотношения (1.4):

Г(^+1) - Г(;„) = -7в |7в| Д*п/2. (2.9)

Соотношение (2.9) применяется обычно при решении задач нестационарного обтекания профиля с угловой кромкой вместо условия (2.7). Получаемые при этом решения оказываются гладкими.

г

Уравнения (2.1), (2.5), (2.6) делают систему уравнений (2.2)-(2.7) нелинейной. Эта система решалась методом итераций, после чего вычислялись интенсивность свободного дискретного вихря Гтп и его комплексная координата по формулам:

Гтп = Г(;п-1) - Г(;га), ^тга(;га) = [¿В(^га) + ¿М+1(;га)]/2.

Комплексные координаты свободных дискретных вихрей Гтт, т = 1, ..., п, для нового момента времени ¿га+1 определялись решением уравнения (1.2):

3. Результаты расчета

Предложенный алгоритм допускает три варианта. Первый (основной) вариант — решение системы уравнений (2.2)-(2.7) для профиля с угловой кромкой ($ = 0) при строгом выполнении условия (2.7) путем задания интенсивности вихревого слоя на последних панелях на профиле квадратичным сплайном вида (2.8). Второй вариант — решение системы уравнений (2.2)-(2.4), (2.6), (2.9). В этом варианте, в отличие от первого, условия (2.5), (2.7) заменены одним условием (2.9), а на последних панелях сохраняется линейная аппроксимация вихревого слоя. Такой подход, как было указано, применяется в большинстве известных алгоритмов решения нелинейных начально-краевых задач обтекания профиля с угловой кромкой. И, наконец, третий вариант — расчет по тем же уравнениям (2.2) - (2.6) с линейной аппроксимацией вихревых слоев на всех панелях без требования выполнения условия (2.7). Этот вариант позволяет строго решать задачу для профиля с острой кромкой ($ = 0). Применение его к расчету профиля с угловой кромкой означает, по существу, замену этой кромки точкой возврата.

Система уравнений (2.2)-(2.6) составляет общую часть всех трех алгоритмов. Эта часть тестировалась с помощью известных решений плоских задач линейной и нелинейной теории тонкого крыла в нестационарном потоке [4, 9].

Следует отметить, что численное решение нелинейных задач имеет свою специфику. Расчет сингулярных интегралов (по соответствующим квадратурным формулам) требует равенства длин элементов на контурах Ь, , что накладывает ограничения на выбор шага по времени. Кроме того, серьезные трудности возникают при расчете нестационарных течений при малых и больших значениях числа Струхаля, а также при больших амплитудах колебаний профиля, вызывающих большие возмущения потока. Все эти особенности учитывались при выборе тестовых задач и основного варианта расчета.

Структура нестационарного течения около задней кромки во многом определяет величину циркуляции скорости вокруг профиля и закон ее изменения с течением времени, а также форму вихревого следа и его интенсивность. Основными локальными характеристиками течения вблизи кромки являются касательные составляющие скорости жидкости на контуре (они равны с обратным знаком величинам 7^, 7(2)), интенсивность сходящих вихрей 7В = к11 7^ + к12 7^ и скорость их схода и>в = (к12 7(2) - к11 7(!))/2. Циркуляция скорости Г является интегральной характеристикой.

Результаты расчета зависимости безразмерных характеристик 7(1);7л2);7в и Г от времени приведена на рис. 2, а ¿Г/^; — на рис. 3. В качестве характерных величин выбирались хорда профиля Ь, скорость потока в бесконечном удалении перед профилем |г^| и период колебаний Т. Расчет проведен для симметричного профиля Кармана — Трефтца,

Рис. 2. Зависимость локальных характеристик течения в угловой кромке и циркуляции скорости от времени (сплошная линия — первый, точки — второй, пунктир — третий варианты расчета).

угловая кромка которого имела угол б = 0.1, а относительная толщина составляла 15%. Профиль совершает поступательные колебания по гармоническому закону с круговой частотой ш и амплитудой 0.01836 от начального момента времени £ = 0 до £ = Т, где Т — период колебаний. Число Струхаля | = 17. Расчет при таких значениях параметров

проводится многими авторами и соответствует условиям эксперимента Брэтта, изучавшего форму вихревого следа за колеблющимся крылом [10].

Точность расчетов контролировалась сходимостью итерационного процесса при решении соответствующих систем уравнений с заданием относительной ошибки вычисления скорости схода вихрей с профиля (величины (1 + ), так как может принимать нулевые значения) на всех шагах по времени. В систематических расчетах относительная погрешность задавалась равной 10-4, что требовало обычно выполнения трех — четырех итераций. Вблизи разрывных режимов течения проводился дополнительный расчет с более высокой точностью, чтобы исключить эффекты погрешности.

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

Сопоставление и анализ результатов расчета позволяют сделать следующие выводы.

1. Для всех трех вариантов алгоритма имеет место гладкая зависимость характеристик течения от времени в начальной стадии движения профиля, что не наблюдается, как правило, в других известных алгоритмах. Эта гладкость — следствие применения метода итераций на каждом шаге по времени, который позволил согласовать начальные данные (выбираемые произвольно) с решением уравнений.

1.5

1.0

<1Г/<И 0-5 о.о

-0.5 -1.0 -1.5

Рис. 3. Зависимость производной циркуляции скорости от времени (слошная линия — первый,

точки — второй, пунктир — третий варианты расчета).

2. Расчет по первому варианту, обеспечивающему строгое выполнение условия (2.7), приводит к существенно негармонической зависимости от времени всех локальных характеристик течения вблизи угловой кромки при гармоническом законе колебаний профиля. Причем решение имеет разрывный характер в моменты времени, в которых происходит смена знака у производной от циркуляции скорости. Такой характер зависимости, как показали дополнительные расчеты, сохраняется для всех значений угла 5 = 0, включая и сколь угодно малые.

3. Расчет тех же характеристик по второму и третьему вариантам алгоритма дает практически гармоническую зависимость от времени, причем результаты расчета по обоим вариантам мало отличаются друг от друга. Исключение составляет лишь скорость схода вихрей с профиля, которая во втором варианте вычисляется по формуле тв = |тв|/2, так как только при этом значении тв соотношение (1.3) переходит в (1.4).

Важно отметить, что решение для полученное по второму варианту алгорит-

ма, не удовлетворяет условию (2.7) в угловой кромке. Таким образом, выполнение условия (2.9) (или его дифференциального аналога (1.4)) не обеспечивает выполнения исходного условия (2.7). Это следует учитывать при формулировке начально-краевой задачи нестационарного обтекания контура с угловой кромкой.

4. Расчет циркуляции скорости, являющейся интегральной характеристикой течения, по трем вариантам алгоритма дает близкие результаты, хотя в первом варианте сохраняется негармонический характер зависимости Г(£) от времени, а производная от циркуляции скорости терпит разрыв в те же моменты, что и локальные характеристики (см. рис. 3).

5. Решение нелинейной начально-краевой задачи нестационарного обтекания профиля с угловой задней кромкой при строгом учете условий в этой кромке в предельном случае бесконечно малого угла 5 не переходит в решение, соответствующее острой кромке (5 = 0).

Что касается формы вихревого следа, то она оказалась мало чувствительной к вариантам алгоритма и хорошо согласуется с экспериментом Брэтта [10] и расчетами по другим алгоритмам [4, 12].

Список литературы

[1] БелоцерковскиЙ С. М., Котовский В. Н., Ништ М. И., Федоров Р. М. Математическое моделирование плоскопараллельного отрывного обтекания тел. Наука,

М., 1988.

[2] БЕТЯЕВ С. К. Эволюция вихревых пелен. В "Динамика сплошной среды со свободными поверхностями". Чебоксары, 1980, 27-38.

[3] Головкин В. А. Нелинейная задача о неустановившемся обтекании произвольного профиля со свободно деформирующимся вихревым следом. Уч. зап. ЦАГИ, 3, №3, 1972, 1-11.

[4] ГОРЕЛОВ Д. Н., КУЛЯЕВ Р. Л. Нелинейная задача о нестационарном обтекании тонкого профиля несжимаемой жидкостью. Изв. АН СССР. МЖГ, №6, 1971, 38-47.

[5] ГОРЕЛОВ Д. Н. Об интегральных уравнениях задачи обтекания профиля. Изв. РАН. МЖГ, №2, 1992, 173-177.

[6] Зовнин А. И. Моделирование начальной стадии развития вихревого следа за телесным профилем. В "Динамика жидкости со свободными границами". Ин-т гидродинамики СО АН СССР, Новосибирск, 60, 1983, 51-59.

[7] Ильичев К. П., Постоловский С.Н. Расчет нестационарного отрывного обтекания тел плоским потоком невязкой жидкости. Изв. АН СССР. МЖГ, №2, 1972, 72-82.

[8] САРПКАЙЯ Т. Вычислительные методы вихрей. Фримановская лекция (1988). Соврем. машиностроение. Сер. А, №10, 1989, 1-60.

[9] Седов Л. И. Плоские задачи гидродинамики и аэродинамики. Наука, М., 1980.

[10] BRATT J. B. Flow patterns in the wake of an ostillating aerofoil. ARC R. and M., No. 2773, 1953.

[11] CHOW C. Y., Huang M. K. The initial lift and drag of animpulsively started aerofoil of finite thickness. J. Fluid Mech., 118, 1982.

[12] Giesing J. P. Nonlinear two-dimensional unsteady potential flow with lift. J. Aircraft, 5, No. 2, 1968, 135-143.

Поступила в редакцию 6 марта 1999 г. в переработанном виде 7 июля 1999 г.