Применение метода конечных элементов для решения задачи обтекания тела вращения под нулевым углом атаки. Характеристика свойства приближения гладких функций интерполяционными многочленами. Расчет лапласиана скалярной функции в декартовых координатах.
Аннотация к работе
Институт проблем механики Российской Академии Наук Потенциальное обтекание тела вращения потоком несжимаемой жидкости. Рассматривается задача о потенциальном обтекании тела вращения потоком несжимаемой жидкости. Приводится численный алгоритм и программы на Фортране. The problem about potential flow of a body of rotation by a flow of incompressible fluid is considered.
Введение
Это шестой препринтт серии объявленной в [1]. Предыдущие препринты этой серии [1-5]. Препринт посвящен описанию програмного комплекса решения задачи об обтекании тела вращения потенциальным потоком несжимаемой жидкости.
I. Численное исследование задачи об обтекании под углом атаки тела вращения потенциальным потоком идеальной несжимаемой жидкости.
Численное исследование устойчивости движения автомобиля при боковом ветре имеет важное значение для безопасности его движения. В настоящее время такие исследования проводятся в аэродинамической трубе обдуванием готового автомобиля или его модели. Натурные эксперименты дороги, кроме того, требуют изготовления опытного экземпляра модели или целого автомобиля. Представляется перспективным заменить натурный эксперимент численным и исследовать аэродинамику автомобиля на стадии проектирования. Для этого требуется выбрать математическую модель и решить численно выписанные уравнения. Поскольку скорости движения автомобиля невысоки по сравнению со скоростью звука и составляют примерно третью его часть, то обычно выбирается модель несжимаемой жидкости. Кроме того, обычно предполагается, что жидкость идеальная и течение является потенциальным. В такой постановке рассматриваемая задача решается панельным методом [6]. Недостатком этого метода
3 является его насыщаемость, т. е. он не использует информацию о гладкости решаемого уравнения Лапласа. В настоящей работе отрабатывается методика построения алгоритма без насыщения на примере тела вращения, обтекаемого потенциальным потоком идеальной несжимаемой жидкости. Построенный алгоритм использует априорную информацию о гладкости решения уравнения Лапласа. Проведенные расчеты показывают перспективность его применения в задаче об обтекании автомобиля потенциальным потоком идеальной несжимаемой жидкости. Количество узлов сетки по сравнению с методом конечных элементов может быть уменьшено в 10 раз при получении в результате такой же точности.
Осесимметричная задача, т.е. обтекание тела вращения под нулевым углом атаки, решалась в [7] методом конечных элементов. Точность полученного решения невысока. В настоящей работе исследуется более общее течение под ненулевыми углами атаки. Предложенный численный алгоритм не имеет насыщения [7] ,т.е. его точность автоматически реагирует на гладкость решения. Это позволяет проводить расчеты с высокой точностью на редкой сетке.
Задача обтекания тела Т потенциальным потоком идеальной несжимаемой жидкости сводится к отысканию потенциала скорости Ф=Ф(х1,х2,х3) - функции гармонической вне обтекаемого тела Т: (1) ?? ? 0, (2) ?? ? 0, ?n
?T
?
(3) grad ? ?V? при x ? ?,
где n - внутренняя нормаль к телу, ДТ - граница тела Т, V? - скорость потока в бесконечности, |х| - длина вектора (х1,х2,х3).
4
Без ограничения общности можно принять, что скорость на бесконечности V? = (1,0,0). В случае шара Т = {х:|х|?1} потенциал имеет вид Ф = х1 х1 /(2|х|3).
?
Как окончательный выход алгоритма интересен не только потенциал скорости, но и давление и вектор скорости, особенно на поверхности обтекаемого тела. Давление вычисляется по формуле, вытекающей из интеграла Бернулли: p??1 ? 2V 2 ? p? ??1 ? 2V 2 , V ? grad?, ?
1 1 где р - давление на бесконечности, ?- плотность, V, V? - модули скорости.
Если в соотношениях (1)-(3) произвести замену Ф=х1 - Ф1 то Ф гармонична в области ?=R3\T: (4) ?n1 ? cos(n,x ), ??
1
(5) ?1 ? 0 при | x |? ?.
Введем криволинейную систему координат (r, ?, ?), связанную с декартовыми координатами (х1,х2,х3) соотношениями: (6) x1=v(r,?)cos?, x2=v(r,?)sin?, x3=u(r, ?).
Если выполняются условия Коши-Римана:
?v 1 ?u ?u 1 ?v ?r r ?? ?r r ??
? ? , ? , то система координат (r, ?, ?) ортогональна, и в этой системе координат Лапласиан скалярной функции имеет вид
Обозначим G область, получаемую меридиональным сечением тела вращения Т. Пусть ?=?(?), ?=u iv, ?=r exp (i?) - конформное отображение круга | ?|?1 на внешность области G, причем центр круга переходит в бесконечность, тогда ? ? r |??(?)|.
Удобно считать, что (r, ?, ?) - сферические координаты, тогда соотношения (6) задают отображение проколотого в центре шара единичного радиуса на внешность рассматриваемого тела Т. В результате замены (6) соотношения (4), (5) принимают вид
(7) ??1 ?|?? | f ( ,?) ? f1( ,?) r?1
?
?
?r
(8) ?1 r?0 ? 0, где f(?,?) получается Из cos(п,х1) заменой (5.6).
Знак в формуле (7) выбран с учетом того, что при отображении (6) внутренняя нормаль к телу переходит во внешнюю.
Произведем еще одну замену: ?2(r, ?, ?)= ?1 (r, ?, ?)- rf1(?, ?), тогда из соотношения (7), (8) получаем
(9) ??2 ? 0, r?1
?r
(10) ?2 r?0 ? 0.
Теперь функция ?2(r, ?, ?) не гармоническая, а удовлетворяет уравнению Пуассона
(11) ??2 ? ??(rf1( ,?) ? F(r,?,?)
? с однородными краевыми условиями (9),(10).
6
Рассмотрим, например, эллипсоид вращения 2 2 2
1 2 3 x x x b2 ? b2 ? a2 ?1? 0.
Тогда:
u(r,?) ? 1?(a ?b)r ? a ?b?cos?, v(r,?) ? 1?(a?b)r ? a ?b?sin?, 2 r
? ?
? ?
2 r
? ?
? ?
?2 ? 1?(a ?b)r ? (a ?b)2 ?2 sin2 ? ? 1?(a ?b)r ? (a ?b)2 ?2 cos2 ?, 4 r 4 r
? ? ? ?
? ? ? ?
? ? ? ? f1( ,?) ? asin? cos?, ?
F(r,?,?) ? a rcos? ?(b ?a)sin2 ? ? 1?(b ?a)r ? a ?b?cos2 ?? ? rsinv cos? .
?
2 2
2 r v
?
? ?
? ?
? ?
? ?
Например, для шара (а=1,b=1) F(r,?,?) ? 2x3 sin? cos? , а решение краевой задачи (9)-(11) есть ? (r,?,?) ? ?r2 ? r?sin? cos? .
2
? ?
? ?
2
? ?
Если ?2 найдено, то для эллипсоида Ф=v(r,?)cos?- ?1(r, ?, ?).
Для дискретизации краевой задачи (9)-(11) применим подход, предложенный в [8]. Таким образом, получаем дискретный Лапласиан в виде h - матрицы: H ? L ? " ?k ?hk , L ? 2l ?1, k?0 l
2 где штрих у знака суммы означает, что слагаемое при k= 0 берется с коэффициентом 1/2, символ ? обозначает кронекерово произведение матриц ?k и hk размера M x M и L x L соответственно. Матрицы ?k размера, M x M, M = тп, получены дискретизацией дифференциального соотношения, зависящего от k:
и рассмотрим интерполяционную формулу g(r) ? ?l? (r)g? , g? ?g(r ), ??1 n
?
(14) l? (r) ? ? " ??p RTP (2r ?1) ? c? r(r ?1)Tm p?0 m m
??p ? m(cos s? ?1) , c? ? p?0 ??p (1? 2p2 ).
?
?
4co
? m?1 m
?
Заметим, что l? (0) ? 0, l? (1) ? 0,т.е. краевые условия (9), (10) 8
? выполнены. Первую, и вторую производные, входящие в соотношение (12), получим дифференцированием интерполяционной формулы (14).
После того, как матрицы ?k, k = 0,1,... ,l, построены для приближенного решения краевой задачи (9)-(11) требуется решить систему линейных уравнений
(15) НФ2=F, где Ф2 и F- векторы, компоненты которых суть значения соответствующих функций в узлах сетки (по ? выбирается сетка из узлов (?k =2?k/L, k=0,1,…,L-1).
В [8] показано, что для обращения h- матрицы справедлива формула
(16) H ?1 ? L ? " ??1k ?hk , L ? 2l ?1, k?0 l
2 т.е. для решения системы линейных уравнений (15) достаточно обратить матрицы ?k, k = 0,1,... ,l, после чего задача сводится к умножению h- матрицы (16) на вектор.
Оценку погрешности построенной дискретизации можно получить стандартными средствами. Отметим только ее качественные особенности. Для дискретизации использовалась интерполяция решения многочленами. Известно [7,9], свойство приближения гладких функций интерполяционными многочленами. Это приближение тем точнее, чем большим условиям гладкости отвечает приближаемая функция. Таким образом, построенный алгоритм не имеет насыщения, т.е. точность полученного приближенного решения тем выше, чем большим условиям гладкости отвечает точное решение.
9
Конкретные расчеты проводились для шара (а=1,b=1) и эллипсоида (а=1,b=0.5) на сетке m=7, n=10, L=9, т.е. состоящей из 630 узлов. Оказалось, что для шара матрица ?0 плохо обращается методом Гаусса. Норма (?0 ?0-1-I) оказалась порядка 10-3. Пришлось применить итерационное уточнение обратной матрицы [10], после чего норма (?0 ?0-1-I) уменьшилась до 10-8. Однако непосредственным умножением матрицы (16) на вектор получить решение не удается. Дело в том, что норка ?0-1 оказалась порядка 1016. Норма остальных матриц ?1, ?2, ?3, ?4, суть 105, 104, 103, 103. Пришлось применить масштабирование с коэффициентом 1010. После этого мантисса полученного приближенного решения совпала с мантиссой точного решения с 7 знаками после запятой (порядок отличался на 10 в соответствии с выбранным масштабным множителем). Для эллипсоида матрица ?0 обращалась хорошо, и итерационное уточнение обратной матрицы не потребовалось.
Нормы матриц ?0, ?1, ?2, ?3, ?4 имели величину порядка 106, 104, 103 ,103, 103, соответственно масштабирование не потребовалось (разумеется, правая часть уравнения Пуассона (11) должна быть вычислена с 6-7 запасными знаками).
Отметим, что для умножения матрицы Н-1 на вектор требуется O(m2 n2LLOGL) операций.
Итак, проведенные расчеты показали, что дискретная задача может быть плохо обусловленной. Поэтому при практическом применении алгоритма следует следить за качеством обращения клеток h - матрицы ?0, ?1, …, ?l. В случае необходимости применять итерационное уточнение обратной матрицы. Также следует следить за нормой полученных матриц ?0-1,?1-1 …,?l-1и применять, если потребуется, масштабирование и вычисление
10 правой части уравнения Пуассона с запасными знаками.
II. Описание программного комплекса.
Решает поставленную задачу для шара программа INT5.
III. Заключение. По поводу получения полных версий описанных программ обращайтесь по электронному адресу: algazinsd@mail.ru или на адрес Института проблем механики РАН, 119526, Москва, проспект Вернадского д.101, к.1.
25
Алгазин Сергей Дмитриевич, Грошев Михаил Владимирович
Численные алгоритмы классической матфизики.
VI. Потенциальное обтекание тела вращения потоком несжимаемой жидкости.
Подписано к печати 2.09.2002. Заказ № 24-2002. Тираж 50 экз. ________________________________________________________
Отпечатано на ризографе Института проблем механики РАН 119526, Москва, пр-т Вернадского, 101