16 ПРИКЛАДНАЯ МЕХАНИКА И ТЕХНИЧЕСКАЯ ФИЗИКА. 2009. Т. 50, N- 6

УДК 532.5.032 : 533.696.5 : 517.54

ОПРЕДЕЛЕНИЕ ФОРМЫ ОСЕСИММЕТРИЧНОГО ТЕЛА,

ОБТЕКАЕМОГО ПОТОКОМ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ,

ПО ЗАДАННОМУ НА ЕГО ПОВЕРХНОСТИ

РАСПРЕДЕЛЕНИЮ ДАВЛЕНИЯ

С. А. Соловьев Научно-исследовательский институт математики и механики им. Н. Г. Чеботарева Казанского государственного университета, 420008 Казань E-mail: [email protected] Разработан способ определения формы осесимметричного тела по заданному распределению коэффициента давления вдоль его меридионального сечения. Учет вязкости осуществлен в рамках модели пограничного слоя. В основу метода положен итерационный процесс, в котором используются решения обратной задачи в плоском случае и прямой задачи для осесимметричного тела. Составлена программа, реализующая итерационный процесс, приведены примеры численных расчетов.

Ключевые слова: обратная краевая задача аэрогидродинамики, вязкая несжимаемая жидкость, пограничный слой, осесимметричное тело, итерационный процесс, панельный метод.

Введение. Обратные краевые задачи аэрогидродинамики (ОКЗА), являющиеся частью общей теории обратных краевых задач, применяются для определения формы крылового профиля по заданному на его поверхности распределению скорости или давления.

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

Одной из основных идеализирующих моделей аэродинамики, упрощающих расчет крыльев летательных аппаратов, является модель плоских сечений, при использовании которой рассматривается не все тело, а лишь его сечения. Такие задачи хорошо изучены как в прямой постановке (см., например, [1–3]), так и в обратной (см., например, [4]).

Разработаны методы, позволяющие получить численно-аналитические решения для различных моделей жидкости: идеальная несжимаемая жидкость (ИНЖ), вязкая сжимаемая жидкость. В случае если применение аналитических методов для расчета профилей невозможно, используются численные методы. Однако при рассмотрении трехмерных тел, имеющих сложную геометрию, модель плоских сечений неприменима. В таких случаях при решении прямых задач используются численные методы с построением сетки, покрывающей всю поверхность тела. При решении пространственных обратных задач также используются численные методы, являющиеся в основном итерационными и поэтому применимые лишь в частных случаях. Например, с помощью модели меридиональных сечений обтекание осесимметричных тел можно рассматривать как двумерную задачу: построение тел Работа выполнена в рамках федеральной целевой программы “Научные и научно-педагогические кадры инновационной России” на 2009–2013 гг.

С. А. Соловьев y Cp x B s r o s 0 L z A V Рис. 1. Постановка задачи:

а — искомое осесимметричное тело; б — распределение коэффициента давления вдоль меридионального сечения Cp (s) (сплошная линия — для осесимметричного тела, штриховая — для симметричного профиля) вращения по заданной хордовой диаграмме распределения скорости [5] или по заданному распределению давления [6], определение формы лопаток турбомашин, расположенных на осесимметричной поверхности тока [7], построение осесимметричного тела по заданному на его поверхности распределению скорости [8].

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

1. Постановка задачи. Непроницаемое осесимметричное тело обтекается осесимметричным потоком ВНЖ с заданной скоростью на бесконечности V (рис. 1,а). Тело задается контуром его меридионального сечения, ось симметрии тела и вектор скорости набегающего потока направлены вдоль оси x. На контуре задано распределение коэффициента давления Cp = Cp (s), s [0, L] (рис. 1,б), дуговая координата s контура меридионального сечения отсчитывается от точки A (s = 0), являющейся точкой разветвления потока, до точки B в задней кромке (s = L). Заданным является также число Рейнольдса Re.

Требуется определить форму меридионального сечения тела.

2. Численно-аналитический способ решения. Эффективным подходом к решению поставленной задачи является подход, предложенный в работе [8], в которой c использованием методов решения обратной задачи для плоского контура и прямой задачи для осесимметричного тела при обтекании в соответствии с моделью ИНЖ построен итерационный процесс решения обратной задачи для осесимметричного тела. Проведенные расчеты показали быструю сходимость процесса (в среднем за 6–8 итераций) с точностью до 106. Используя методы решения обратной плоской задачи и прямой осесимметричной задачи по модели ВНЖ, построим аналогичный итерационный процесс.

(0) 1. Выберем начальное приближение распределения коэффициента давления Cp1 (s) для (0) решения ОКЗА в плоском случае: Cp1 (s) = Cp (s), k = 0.

(k) (k) 2. По распределению Cp1 (s) найдем контур Lz, решив ОКЗА для симметричного профиля.

18 ПРИКЛАДНАЯ МЕХАНИКА И ТЕХНИЧЕСКАЯ ФИЗИКА. 2009. Т. 50, N- (k) 3. В случае если построенный контур Lz оказался разомкнутым, замкнем его с помощью метода квазирешения [4]. При этом если контур имеет самопересечения, то форма профиля подправляется до однолистного.

4. Определим Cp2 (s) как решение прямой краевой задачи для осесимметричного теk) ла, меридиональным сечением которого является контур Lz с заданной скоростью V набегающего потока.

5. Вычислим невязку (k) (s) = Cp (s) Cp2 (s).

6. Следующее приближение распределения коэффициента давления для плоской ОКЗА определим по формуле Cp1 (s) = Cp1 (s) + (k) (s), где [0, 1] — коэффициент релаксации.

В качестве критерия завершения итерационного процесса используется условие (k) (s) < µ (µ — заданное положительное число). В этом случае полученное решение с заданной точностью совпадет с исходным распределением коэффициента давления.

Так как в общем случае обратная задача может не иметь решения в классе замкнутых контуров, то вторым условием завершения процесса является следующее условие:

Cp2 (s) Cp2 (s) < µ. В этом случае распределение коэффициента давления на полученном профиле может отличаться от заданного.

Следует отметить, что при использовании метода квазирешений для решения плоской ОКЗА учитывалось только условие замкнутости контура профиля, при этом скорость набегающего потока, полученная в процессе решения, могла не совпадать с заданной скоростью V. Прямая краевая задача аэрогидродинамики решалась для изначально заданной скорости набегающего потока. Если полученный контур имеет самопересечения, то однолистности профиля можно добиться, увеличивая угол наклона касательной к кривой зависимости Cp1 (s) в окрестности задней кромки (штриховая линия на рис. 1,б).

3. Решение обратной краевой задачи для симметричного профиля. При решении обратной краевой задачи в плоском случае требуется определить форму симметричного крылового профиля Lz (рис. 2,а), безотрывно обтекаемого потоком ВНЖ, по заданному на его поверхности распределению коэффициента давления Cp = Cp (), [0, 2L] (рис. 2,б). Дуговая координата контура профиля отсчитывается от точки B, так чтобы при возрастании область течения оставалась слева. Точка A с дуговой абсциссой = = L является точкой разветвления потока.

При решении ОКЗА по модели ПС задача сводится к нахождению полутела вытеснения CB AB C по заданному на участке B AB распределению скорости V = V () [4].

Рис. 2. К решению ОКЗА в плоском случае:

а — физическая плоскость z; б — распределение коэффициента давления по контуру С. А. Соловьев Зависимость V () определяется по заданной функции Cp () в формуле, следующей из интеграла Бернулли:

Задача отыскания полутела вытеснения сходна с основной ОКЗА для профиля в потоке ИНЖ. В случае если искомый профиль и обтекающий его поток симметричны, метод решения обратной задачи упрощается.

Так как течение является безвихревым, то существует потенциал скорости (x, y) и связанная с ним функция тока (x, y). Введем в рассмотрение комплексный потенциал течения w(z) = (x, y) + i(x, y). Из предположения непроницаемости контура полутела вытеснения следует, что на нем = 0. Тогда, полагая ( ) = 0 на поверхности полутела, получим Обозначим 1 = (2L). Так как рассматривается симметричное обтекание симметричного тела, то циркуляция = 0.

Введем в рассмотрение каноническую плоскость = r ei, в которой области обтекания полутела вытеснения соответствует внешность круга единичного радиуса || 1с разрезом от точки = 1 до бесконечности. Для взаимно однозначного конформного отображения этих областей потребуем соответствия бесконечно удаленных точек плоскостей z и, а также соответствия точки = 1 точке z = 0. При этом граница единичного круга || = 1 перейдет в границу полутела, а верхний и нижний берега разреза — в линии B C и B C соответственно. Комплексный потенциал симметричного потока в канонической плоскости имеет вид где u0 — модуль скорости на бесконечности в плоскости ; C0 — постоянная.

Приняв = ei и учитывая, что потенциал скорости () = Re [w(ei )] в точках = и = 1 известен: (0) = 1, () = 0, получим систему уравнений для нахождения неизвестных параметров, решив которую найдем C0 = 2u0 = 1 /2. В итоге на границе L имеем Так как в плоскостях z и значения потенциалов скорости () и () в соответствующих точках совпадают, то из формул (1) и (2) следует зависимость = (), [0, 2].

Введем в рассмотрение вспомогательную, аналитическую в области || 1 функцию и получим граничное значение ее действительной части () = ln |ds/d| при = ei.

Так как полутело вытеснения не должно иметь изломов в точках B и B, т. е. внутренний по отношению к области течения угол равен, а в окрестности образа этих точек в плоскости контур L имеет угол /2, то в этой окрестности z() (1 1/)2. Следовательно, функция () имеет логарифмическую особенность в точке = 1. Устранив эту особенность, получим функцию 20 ПРИКЛАДНАЯ МЕХАНИКА И ТЕХНИЧЕСКАЯ ФИЗИКА. 2009. Т. 50, N- Зная зависимость = (), найдем P () = Re [()] на границе L и по формуле Шварца восстановим функцию () с плотностью P ():

Функция z() определяется по формуле Подставив в (3) предельное значение = ei, получим параметрические уравнения искомого контура полутела вытеснения Для того чтобы определить форму самого профиля, надо сместиться по нормали к границе полутела вытеснения на участке B AB внутрь этого полутела на расстояние, равное толщине слоя вытеснения 1 (), которая может быть найдена путем расчета ПС любым известным методом (см., например, [2, 9–11]).

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

Существование физически реализуемого профиля подразумевает замкнутость его контура. Для модели ПС замкнутость контура профиля означает разомкнутость полутела вытеснения на величину 1 = 21 (0) = 21 (2L). Контур крылового профиля будет замкнут, если имеет место равенство Для того чтобы условия разрешимости были выполнены, применим метод квазирешения [4]. Кроме того, полученный контур не должен иметь самопересечений.

4. Решение прямой краевой задачи для осесимметричного тела. При решении прямой краевой задачи требуется определить распределение коэффициента давления вдоль контура меридионального сечения Cp (s), s [0, L] для осесимметричного тела, обтекаемого осесимметричным потоком ВНЖ с заданной скоростью на бесконечности V.

Ось симметрии тела и вектор скорости набегающего потока направлены вдоль оси x (см.

рис. 1,а).

Решение прямой задачи в осесимметричном случае проводится панельным методом (см., например, [1, 12]). В качестве панелей выбираются кольца, расположенные вдоль тела вращения. Перейдем к цилиндрическим координатам x,,, связанным с декартовыми x, y, z соотношениями y = sin, z = cos, и рассмотрим меридиональное сечение при = 0. Тогда на меридиональном сечении тела панелями являются последовательно соединенные отрезки {xj, j }, j = 1, N + 1, где N — число панелей. В качестве контрольных точек выберем середины панелей с координатами {xjc = (xj +xj+1 )/2, jc = (j +j+1 )/2}, С. А. Соловьев Рис. 3. Форма меридионального сечения осесимметричного тела, полученная при решении прямой задачи:

1 — контур полутела; 2 — контур осесимметричного тела Полагается, что на поверхности тела распределены гидродинамические особенности:

источники и стоки. Обозначим через qj интенсивность распределенных источников j-й панели (qj > 0 — источник, qj < 0 — сток), интенсивность на каждой панели будем полагать постоянной вдоль ее длины. Значения qj определяются в ходе решения задачи.

Введем следующие обозначения: j — угол между j-й панелью и положительным направлением оси x (aj = tg j ), rkj — расстояние от текущей точки интегрирования j-й панели с особенностями до контрольной точки k-й панели (рис. 3).

Потенциал скорости (x, ) складывается из потенциала невозмущенного потока = V x и суммы потенциалов j (x, ), индуцируемых гидродинамическими особенностями на панелях. Потенциал в k-й точке принимает вид Выражение (4) записано в общем виде и может быть использовано при любом выборе панелей и способе распределения особенностей. Нормальная и касательная составляющие скорости выражаются через производные от потенциала:

Здесь nxk = sk = sin k ; nk = sxk = cos k.

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

Тогда выражение для нормальной составляющей скорости принимает вид где K(mkj ) — полный эллиптический интеграл первого рода; E(mkj ) — полный эллиптический интеграл второго рода; mkj = 4k j / (k j )2 + (xk xj )2. Таким образом, интеграл в соотношении (7) имеет особенность в случае, когда точка k совпадает с точкой j. Для того чтобы избежать сингулярности интеграла при k = j, контрольные точки необходимо сместить в поток на расстояние от контура профиля. Координаты этих точек равны xk = xk + nxk, k = k + nk [3]. Сами интегралы находятся численно с использованием квадратурных формул. Опыт расчета обтекания сферы потоком ИНЖ показал, что совпадение с аналитическим решением ( < 108 ) достигается на расстояниях от контура, составляющих 0,01 1,00 длины панели. Во всех расчетах значение равно 0, длины панели.

При использовании панельного метода в случае ИНЖ значения qj вычисляются из условия непротекания потока через поверхность тела При обтекании вязкой жидкостью в соответствии с моделью ПС используется условие непротекания через поверхность полутела вытеснения. В этом случае нормальная составляющая скорости на поверхности тела не равна нулю:

Значение Vn можно определить из уравнения неразрывности для пограничного слоя осесимметричного тела [12, 13] где Ue (s) — распределение скорости на поверхности полутела вытеснения при его обтекании потоком ИНЖ; 1 (s) — распределение толщины вытеснения ПС; R(s) — радиус меридионального сечения.

Соотношение (8) представляет собой систему линейных алгебраических уравнений относительно неизвестных qj, используемую для нахождения распределения скорости на поверхности тела. Искомое распределение скорости на поверхности полутела вытеснения Ue (s), а также зависящее от него распределение 1 (s) входят в правую часть (8), поэтому для определения Ue (s) строится следующий итерационный процесс.

1. Начальное приближение Ue (s) (n = 0) берется из решения задачи в рамках модели ИНЖ.

2. С учетом (9) решается система (8).

3. После определения потенциала (n) (xk, k ) по формуле (6) находятся компоненты скорости Vs (xk, k ) вдоль контура тела.

5. Принимается n = n + 1 и выполняется переход на шаг 2.

С. А. Соловьев Процесс продолжается до тех пор, пока не будет выполнено условие |Ue (s) Ue (s)| < µ. Тогда Ue (s) является искомым распределением скорости на поверхности полутела вытеснения. Распределение коэффициента давления на поверхности тела можно найти по формуле, следующей из интеграла Бернулли:

При расчете обтекания для осесимметричного ПС законы сохранения массы и количества движения записываются в виде где u, v — компоненты вектора скорости; U — скорость на внешней границе пограничного слоя; — кинематическая вязкость; r — расстояние от оси симметрии до точки в пограничном слое. При k = 0 система уравнений (10) описывает плоский ПС, при k = 1 — осесимметричный ПС. Задаются следующие граничные условия:

Поскольку в системе (10) пограничный слой является тонким, значения r(s, n) можно заменить на R(s), что позволяет существенно упростить дальнейшие вычисления. Такая замена может ограничить область применимости уравнений ПС (например, в задней кромке, где толщина пограничного слоя достаточно велика, а радиус меридионального сечения стремится к нулю). Однако численные расчеты, проведенные в данной работе, а также расчеты и результаты сравнения с экспериментальными данными в работах [12, 13] показали, что это допущение можно использовать и при очень малых значениях R(s), а значение в последних точках получать экстраполированием.

Уравнения осесимметричного ПС могут быть сведены к уравнениям плоского ПС с помощью преобразования Манглера — Степанова (см., например, [9, 10]):

В дальнейшем расчет ПС проводится для преобразованных координат интегральным методом Кочина — Лойцянского (см., например, [2, 9, 10]) или Эпплера (см., например, [11, 12]).

Коэффициент сопротивления вычисляется с помощью аналога формулы Сквайра — Юнга для осесимметричного тела [12, 13]:

(2 (s) — толщина слоя потери импульса; V — объем тела; U — скорость на бесконечности; H12 = 1 /2 — формпараметр). В отличие от плоского случая в качестве коэффициента сопротивления осесимметричного тела принимается величина max Cx (s), так как в задней кромке R(s) 0.

5. Тестовые расчеты. Для тестирования изложенного способа проведены численные расчеты, результаты которых представлены на рис. 4, 5. В качестве исходных были взяты распределения коэффициентов давления, полученные при решении прямой задачи для заданных тел (сплошные линии на рис. 4). Во всех расчетах Re = 106. Все рассмотренные 0, Рис. 4. Распределение коэффициента давления вдоль меридионального сечения тела, полученное в результате тестовых расчетов:

а — тело толщиной менее 15 %, б — тело толщиной более 15 %; сплошные линии — исходное распределение, штриховые — распределение, полученное в результате решения ОКЗА, точки — решение задачи в CFD-пакете “Fluent” 0, Рис. 5. Форма меридионального сечения тела, полученная в результате тестовых а — тело толщиной менее 15 %, б — тело толщиной более 15 %; сплошные линии — исходная форма, штриховые — форма, полученная в результате решения ОКЗА тела имеют острую заднюю кромку с внешним углом, равным 2, и обтекаются безотрывно. Для тел, имеющих толщину менее 15 % длины хорды (рис. 5,а), итерационный процесс сходится быстро (в среднем за 4–6 итераций) с точностью до 106. Полученное в результате решения меридиональное сечение и соответствующее ему распределение коэффициента давления показаны на рис. 4, 5 штриховыми линиями.

Наряду с тонкими телами рассмотрены тела большей толщины (сплошная линия на рис. 5,б). Для безотрывно обтекаемого тела в результате решения прямой краевой задачи аэрогидроупругости получено распределение коэффициента давления (сплошная линия на рис. 4,б). Однако уже на третьем итерационном шаге при решении ОКЗА для симметричного профиля расчет ПС показал, что обтекание будет отрывным. Отрыв происходит в хвостовой части тела, что обусловлено наличием отрицательного градиента скорости.

Следовательно, получить меридиональное сечение исходного тела невозможно, так как ОКЗА с учетом вязкости имеет решение только в предположении безотрывности обтекания. Одним из способов преодоления этой трудности является решение плоской ОКЗА с использованием модели ИНЖ. Иными словами, по заданному распределению коэффициента давления строится не полутело вытеснения, а симметричный профиль, обтекаемый потоком ИНЖ. Полученное в итоге меридиональное сечение и соответствующее ему распределение коэффициента давления показаны на рис. 4,б, 5,б штриховыми линиями. При С. А. Соловьев 0, Рис. 6. Результаты решения задачи определения формы осесимметричного тела:

а — распределение коэффициента давления вдоль меридионального сечения тела (сплошная линия — заданное распределение, штриховая — распределение, полученное в результате решения ОКЗА); б — форма меридионального сечения тела этом итерационный процесс сходится за большее число итераций (20 итераций). Это вызвано тем, что влияние вязкости наиболее существенно в хвостовой части тела. Решение с учетом наличия ПС позволяет более точно определить угол раствора и толщину задней кромки. Применение модели ИНЖ, для того чтобы решение ОКЗА удовлетворяло условию замкнутости, приводит к снижению точности этого решения, в то время как учет вязкости при решении обратной плоской задачи обусловливает более быструю сходимость итерационного процесса с большей точностью.

Для оценки достоверности полученных результатов проведены расчеты в CFD-пакете “Fluent”. Для этого в сеточном генераторе “Gambit” построены регулярные прямоугольные сетки. Расчеты проводились с использованием модели турбулентности Спаларта — Алмараса. Распределение коэффициента давления, полученное для примера, представленного на рис. 5,б, показано на рис. 4,б точками. Из анализа результатов сравнения расчетов следует, что распределения коэффициентов давления хорошо согласуются ( 108 ).

Для исходного тела коэффициент сопротивления Cx = 0,026, для спроектированного тела Cx = 0,025. Близость значений коэффициентов сопротивления, полученных в проведенных расчетах, объясняется тем, что в формуле (11) в качестве характерного размера принята площадь S = V 2/3. Иными словами, при близкой геометрии тел изменение аэродинамических характеристик R, 2, Ue в (11) компенсируется изменением площади меридионального сечения.

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

На рис. 6 сплошной линией показано исходное распределение коэффициента давления;

контур меридионального сечения и соответствующее ему распределение коэффициента давления, полученные в результате решения ОКЗА, представлены штрихпунктирной и штриховой линиями соответственно. Коэффициент сопротивления Cx = 0,026. Отличие значения коэффициента давления в хвостовой части от исходного значения обусловлено выполнением условий разрешимости. Так как в общем случае обратные задачи могут не иметь решения, то можно утверждать, что в этом случае для заданного распределения коэффициента давления и заданной скорости набегающего потока не существует физически реализуемого тела.

Заключение. Таким образом, решение задачи построения осесимметричного тела, обтекаемого потоком ИНЖ, обобщено на случай учета вязкости по модели ПС, что позволяет определить значение Cx. Анализ результатов проведенных численных расчетов показал, что разработанный итерационный способ решения ОКЗА для осесимметричных тел, обтекаемых потоком ВНЖ, является эффективным. Вследствие использования в процессе решения метода решения обратных задач и численного метода решения прямых задач аналитическая оценка сходимости представляет собой сложную задачу. Однако проведенные численные расчеты позволяют сделать вывод, что предложенный итерационный метод сходится, причем достаточно быстро (в среднем не более чем за 20 итераций). Применение панельного метода для моделирования вязкого обтекания при решении прямой задачи показало его эффективность, что обеспечивает высокую скорость работы программ, реализующих разработанный способ.

Автор выражает благодарность Н. Б. Ильинскому за предложенную тему исследования и полезные советы.

ЛИТЕРАТУРА

1. Флетчер К. Вычислительные методы в динамике жидкостей: В 2 т. Т. 2. Методы расчета различных течений / Под ред. Л. И. Турчака. М.: Мир, 1991.

2. Лойцянский Л. Г. Механика жидкости и газа. М.: Наука. Гл. ред. физ.-мат. лит., 1987.

3. Бенерджи П. Метод граничных элементов в прикладных науках / П. Бенерджи, Р. Баттерфилд. М.: Мир, 1984.

4. Елизаров А. М. Обратные краевые задачи аэрогидродинамики / А. М. Елизаров, Н. Б. Ильинский, А. В. Поташев. М.: Наука, 1994.

5. Киселев О. М. Построение тела вращения по заданному на нем распределению скорости // Изв. вузов. Авиац. техника. 1959. № 2. С. 20–24.

6. Этерман И. И. Определение поверхности тела вращения по заданному распределению давления // Докл. АН СССР. 1947. Т. 56, № 4. С. 351–353.

7. Вячкилев О. А., Ильинский Н. Б., Исмагилова Г. Р. и др. Обратная краевая задача для решетки профилей, расположенной на осесимметричной поверхности тока в слое переменной толщины // Журн. вычисл. математики и мат. физики. 1996. Т. 36, № 11. C. 134–142.

8. Ильинский Н. Б., Марданов Р. Ф., Соловьев С. А. Комбинированный метод решения обратной краевой задачи аэрогидродинамики для осесимметричного тела // Журн. вычисл.

математики и мат. физики. 2008. Т. 48, № 7. С. 1294–1308.

9. Лойцянский Л. Г. Ламинарный пограничный слой. М.: Физматгиз, 1962.

10. Шлихтинг Г. Теория пограничного слоя. М.: Наука, 1974.

11. Eppler R. Airfoil design and data. Berlin: Springer-Verlag, 1990.

12. Lutz Th. Berechnung und optimierung subsonisch umstromter prole und rotationskorper.

Dusseldorf: VDI Verlag, 2000.

13. Hess J. L. On the problem of shaping an axisymmetric body to obtain low drag at large Reynolds number // J. Ship Res. 1976. V. 20. P. 51–60.