Бабенко - Основы численного анализа (947491), страница 148
Текст из файла (страница 148)
На сУГ введем локальные координаты (х, е). Поскольку хя = х, х = е(х) соед, Заметим, что на дТ никакие условия на ф и яе не накладываются (см. п. 4). Изложим решение этой задачи методом конечных элементов. Можно сразу же сказать, что мы пе получим эффективного алгоритма, так как не сможем воспользоваться гармоничностью функции зз в полной мере. Однако, чтобы перейти к более трудной нелинейной задаче об обтекании тела я:жнмаемой жидкостью, необходимо предварительно решить задачу о течении несжимаемой жидкости, причем таким методом, который без труда переносится на более общий случай, Именно поэтому выбран вариационный метод конечных элементов. Поскольку рассматриваемая задача являетгя исгледовательской, метод решения системы линейных уравнений выбирался из соображений простоты реализации, а не максимальной эффективности. Описываемые ниже вычисления проделаны Л.
И. Аптекаревым и Л, Р. Волевичем (см, также (223-2261). Ограничимся случаем осесимметри.п|ых течений. Учитывая, что И = (1, О, О), течение будем считать осесимметричным., если тело является телом вращения относительно оси хя. В цилиндрических координазах (х,, у, д), хя = х, хз = усове, хз = уе!и О потеппиал сз не зависит от О, и поэтому вариационную задачу нужно рассматривать в классе функций, инвариантных относительно вращении вокруг оси хя, т, е, функций еь зависящих только от х, у. Пусть дуга Г, от вращения которой вокруг оси хя получается поверхность тела ОТ, является графиком функции у — — е(х). Ооозначим через Те область, ограниченную дуяой Г и отрезком оси х, на который она опирается, и пусть По .=- (х, у: — оо < х < оо, О < у < оо) 'з То Поскольку в цилиндрических 724 Глава 10.
Некотяорые вопросы часленного решснил краевых задач хз = з(х) сбпО на дТ, то коэффициенты первой квадратичной формы поверх- ности запишутся в виде Š—. ~ ( — ~) —.- 1-~- зв, Е = ~( — '~) — —. е (х), С =. ~( — — ) = О. с=1 т=т Поэтому мера Лебега на дТ примет вид с1о = чтЕà — Сгдхт1О = в(1 —: + з ) т1хдО. Вектор нормали найдется по обычным правилам; тт = хт х тг тт'2 хго/(г . х го„где г — радиус-вектор точки поверхности. Если произвести выкладку и учесть направление нормали, то получим тг = (е'(х), — совО, — з!пО)/(1, е' )" Таким образом, считая, что дуга Г опирается на интервал ( — 1, 1), найдем 1 Ьь~<1о =- 2тт з( г (х)з(х) ~И(хт в(х))с(х. сп Итак, вариационное уравнение (54) запишется я виде 1 Л вЂ” — — — — 1 д дх т1з1 — 1 и (х) з(х) ф(х., з(х)) дх = О, тдтр дт(т дтр дйз Ьхд дрдр) / (55) пс где т)т — произвольная функция, для которой (тр й) т ~ В рстхс17О ч 00, При численном решении задачи заменим область йо на ограниченную подобласть йя, содержащую круг достаточно большого радиуса и ограниченнуто дугами 1 и, Г и двумя отрезками оси и = О (рис.
5). Учитывая, что го = — 0((х~ г), на Гя поставим нулевое граничное условие. Итак, если тр — искомая гармоническая функция, то она будет удовлетворять уравнению (55), в котором йе нужно заменить -1 Рис. 5. 0бласть, в котоРой па йи и считать, что В б Иг(йв) и Ет~г решается задача обтекания =.. О. Кроме того, тх( = О, а компонента гра- та ницы (АВ) Г С,'СГУ) свободна от постановки граничных условий. Приближенное решение этой задачи будем трактовать как приближенное решение исходной задачи.
Перейдем к неформальному описанию численного алгоритма. Первын шаг, который надлежит сделать при чпштенной реализации метода конечных элементов, это произвести триангуляцию области. Для этого нужно в области ввести сетку и после соединения узлов сетки отрезками полу.чить триангуляцию. Никаких универсальных методов построения сеток нет. Вопрос этот 725 Ь" 2. Вориациоиные мееподы решет л краевых задач трудный, и на практике построение хороших сеток достигается с помощью многочисленных предварительных расчетов.
Обычно поступают следующим образом: строят гомеоморфизм области на некоторую каноническую область. В двумерном случае для односвязной области в кжестве канонической области принимают прямоугачьник. На канонической области строят регулярную сетку и затем с помощью гомооморфизма переносят ее на исходную область. Требуемый гомеоморфизм часто строят как конформное или квазиконформное отображение области на прямоугольник (27). Если отображается область с гладкой границей., то при таком подходе возможны значительные искажения при переходе от образа (прямоугольника) к прообразу (области) Построив сетку и произведя триангуляцию, необходимо занумеровать вершины.
Проце10ра нумерации неформальная и очень важная, поскольку способ ну.мерашеи определяет положение ненулевых элементов в матрице системы (32). В общей ситуации процедура нумерации нетривиальна. В рассматриваемой задаче указанные вопросы решанзтся весьма просто, если счнтатеч что тело является эллипсоидом вращения вокруг оси хь Ниже это и будет предполагаться. Таким образом. ,р:х +уХЬ =1,д>0).
Введя эллиптические координаты при Ь < 1: х = ссйб гово, р = св!е6эшг1, где с = 11 сйбо, бо = агсеЬЬ, и приняв, что Гн — дуга эллипса 6 = В, получим гомеоморфизм области Нл на прямоугольник 1л .= (Д, 6: Хо < 6 < В., 0 ( 0 ( г). В 1л введем сетку (Ч, т1,), полагая Е=чео-~а~ ~, ой =я —, (( —,'(Хх -1), г=1, 2,..., ЛХШ, (56) где (6) — дробная часть О, причем мы принимаем (1) 1. Константа а в результате экспериментов была принята равной 0,12. Так м образом, Л =.. Ео — ' -~- а(АХ вЂ” 1), и при Ь = 172 полуоси эллипса, на котором лежит дуга Гл, будут равны 80,83 и 80,82 соответственно.
Если тело -- шар, то место эллиптических координат брались полярные координаты (г, и) и прин малосли что г, = о' (1=1,2,...,61,о>1. Таким образом, в плоскости (х, и) координаты узлов сетки Р, (х„р.) бучут равны х; = сс11Е сов ц„ш — — св116, ыпгЬ (1 =. Н 2, ..., ЛХУ), и формулы (56) устанавливают их нумерацию. В случае сферических координат способ нумерацепе указан на рис. 6, и он же остается таковым для Эллиптических координат, толькО мсняется характер координатных линий. Заметим, что все узлы распадшотсч на М слоев по а Р точек в каждом етое, причем стой, для которого 1 (1 < а, .лежит на Г.
рис. 6 Н .мерация узлов сетки Чтобы получить триангуляцию области, произведем ее в параметрической пл~юкости -- в прямоугольнике 1л, а затем перенесем на плоскость (х, й), но узлы будем соединять прямолинейнылп~ отрезками. Одна из возможных тришегуляпий (называемая лля удобства правой) описшеа перед формулировкой зада еи 6.
Можно использовать иную триангуляцию (левую), проведя диагональ Ро еР,+ы в четырехугольнике РоР,+г1ххем+~Ргеы (см. рис. 6) общего 726 Глава 10. Неко>порыв оопросы численного решения краевых вас)ач положения. Потенциал вы является нечетной функцией х, и для того, чтобы сохранить это свойство, а численном реснении триангуляцию целесообразно делать так, чтобы картина была симметричной относительно оси р. т.е.
изменить триангуляцию с левой на правую на некоторых лучах. В каждом из треугольников триангуляции делается линейная интерполяция, и поэтому базисная функция е,(х, р), отвечалнцая узлу Р,, будет определяться по формулам задачи б, если триангуляция правая; при левой триангуляции нужно сделать очевидные изменения и конструкции носителя базисной функции. Если же симплекс примыкает к дуге Г либо оси р =- О, то базисная функция, отиечающая узлам., лежащим на Г либо на отрезках оси р = О, эидоизменяется, и соответствующие изменения может внести сам читатель. Вычисление интегралов (е>о ес) выполняется элементарно и в явном виде; надлежащие формулы имев ются в задаче 7.
Исключение составляют треугольники, примыкающие к дуге Г, и пусть сЛЯЯЯг — один из таких треугольников (рис. 7). Из него удаляется защтрихованный криволинейный сегмент Е. Если следовать предписаниям теории, интеграл нужно брать по кривозинейно>су треугольнику >Л(Х>Г>Яв с, Е, а можно >'>с ВЫЧИСЛятЬ НитЕГраЛ И ПО СЗХ>>>1Х>Х>>В. Прн ПрИНятОй ИН- терполяции (см. об этом ни>ко), как показали численные эксперименты, уточнения ращения не про>>ссадит, если вычислять интсгралы по кривОлинейному треугальнику. Вычислим правые части системы линейных урааненнй. Из формулы (55) заключаем, что если узел Р> Г, то в уравнении, отвечающем функции е,(х, р), правая часть равна Рис.
7. Криволиней пый граничный эле мент > о — -- / в'(х)в(х)е>(х, в(х))с)х, Х =- 1, 2, ..., >У. — 1 (о7) Если Х ( > < 2дс, то вклад от этих узлов в правую часть будет равен нулю при замене кривой Г на ломаную Р,... Рн. В противном случае вклад будет третьего порядка малости, и мы им пренебрежем. При 7' ) 2Л' а> = О. Интеграл (57) вычисляется элементарно, поскольку после интегрирования по частям получим 1 Х, А а — —. — — / в (х) — е (х, в(х))с)х, 2,/ с)х — 1 а этот интеграл вычисляется я яаном виде.
Таким образом, система (32) и дан- ном случае приобретает вид А>ь»> =а>, с=.. 1>2,...> ЛХЖ, (58) с=> где >сь — приближенное значениепотенциала р а узле Рь. срормсбаы для коэффициентов системы можно выписать самостоятелып>. Запишем систему (58) в иной форме. Введем векторы Л> = (Гыс» О ..., тн» О н)' О = 1, 2, .... ЛХ), с>, = (ав>О .>1+>, ..., асей О, г;)' (Х .—...
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
