авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 | 3 | 4 |
-- [ Страница 1 ] --

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РФ

КЕМЕРОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

КАФЕДРА ЮНЕСКО ПО НОВЫМ ИНФОРМАЦИОННЫМ ТЕХНОЛОГИЯМ

К.Е. Афанасьев, А.М. Гудов

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

В ЧИСЛЕННЫХ РАСЧЕТАХ

Учебное пособие

Кемерово 2001

ББК B253.31я73+B192.1я73+3811я73

А94

УДК 532.5.031+519.6+681.3

Печатается по решению редакционно-издательского

Совета Кемеровского университета

Рецензенты: доктор физ.-мат. наук, профессор В.П. Житников, кафедра вычислительной техники и информацион ных технологий Кузбасского государственного тех нического университета А94 Афанасьев К.Е., Гудов А.М. Информационные технологии в численных расчетах: Учебное пособие. – Кемерово: КемГУ, 2001.- 204с.

ISBN 5-8353-0063- В пособии рассмотрены основные вопросы гидродинамики идеальной жидкости со свободными границами. Дано описание численных методов граничных элементов для ре шения таких задач.

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

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

А Без объявл ББК ЛР№ © К.Е. Афанасьев А.М. Гудов, ISBN 5 -8353-0063-8 © Кемеровский государственный университет, Оглавление ОГЛАВЛЕНИЕ ВВЕДЕНИЕ Г л а в а 1. НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ 1.1. НАЧАЛЬНЫЕ СВЕДЕНИЯ О КРАЕВЫХ ЗАДАЧАХ.............................................. 1.2. ОСНОВНЫЕ УРАВНЕНИЯ ТЕОРИИ ИДЕАЛЬНОЙ ЖИДКОСТИ.......................... УРАВНЕНИЕ НЕРАЗРЫВНОСТИ УРАВНЕНИЕ ДВИЖЕНИЯ ИНТЕГРАЛЫ ДВИЖЕНИЯ Г л а в а 2. ПОСТАНОВКА ЗАДАЧ СО СВОБОДНЫМИ ГРАНИЦАМИ, РЕШАЕМЫХ ПРИ ПОМОЩИ ПАКЕТА “AKORD” 2.1. ВЕЩЕСТВЕННЫЕ ПЕРЕМЕННЫЕ................................................................... ПЛОСКИЙ СЛУЧАЙ ОСЕСИММЕТРИЧНЫЙ СЛУЧАЙ 2.2. КОМПЛЕКСНЫЕ ПЕРЕМЕННЫЕ. ПЛОСКИЙ СЛУЧАЙ................................... 2.3. ПРИВЕДЕНИЕ КРАЕВЫХ ЗАДАЧ К БЕЗРАЗМЕРНОМУ ВИДУ........................... Г л а в а 3. МАТЕМАТИЧЕСКИЕ МЕТОДЫ И АЛГОРИТМЫ, ИСПОЛЬЗУЕМЫЕ В ППП «AKORD» 3.1. КМГЭ НА ОСНОВЕ ИНТЕГРАЛЬНОЙ ФОРМУЛЫ КОШИ.............................. 3.2. МОДИФИЦИРОВАННЫЙ КМГЭ.................................................................... 3.3. МГЭ НА ОСНОВЕ ТРЕТЬЕЙ ФОРМУЛЫ ГРИНА............................................ ПЛОСКИЙ СЛУЧАЙ ОСЕСИММЕТРИЧНЫЙ СЛУЧАЙ 3.

4. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ 3.5. АЛГОРИТМЫ ПОСТРОЕНИЯ СВОБОДНЫХ ГРАНИЦ....................................... СТАЦИОНАРНАЯ ЗАДАЧА НЕСТАЦИОНАРНАЯ ЗАДАЧА 3.6. РАБОТА С ГРАФИЧЕСКИМИ ОБЪЕКТАМИ..................................................... БАЗОВЫЕ ГРАФИЧЕСКИЕ ОБЪЕКТЫ НА ПЛОСКОСТИ БАЗОВЫЕ ПРОСТРАНСТВЕННЫЕ ГРАФИЧЕСКИЕ ОБЪЕКТЫ 3.7. ПОСТРОЕНИЕ СЕТОК В МЕТОДАХ ГРАНИЧНЫХ ЭЛЕМЕНТОВ...................... ПОСТРОЕНИЕ ДВУМЕРНОЙ ЧЕТЫРЕХУГОЛЬНОЙ СЕТКИ ПОСТРОЕНИЕ ПРОСТРАНСТВЕННОЙ ТРЕУГОЛЬНОЙ СЕТКИ Г л а в а 4. ПАКЕТ ПРИКЛАДНЫХ ПРОГРАММ “AKORD” 4.1. ПРЕПРОЦЕССОР.............................................................................................. 4.2. РЕШАТЕЛИ...................................................................................................... 4.3. ПОСТПРОЦЕССОР........................................................................................... Г л а в а 5. ОТ ОТДЕЛЬНЫХ ПРОГРАММ К ИНТЕГРИРОВАННОЙ ИНФОРМАЦИОННОЙ СИСТЕМЕ ПОДДЕРЖКИ ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА 5.1. СТРУКТУРА ИНФОРМАЦИОННОЙ СИСТЕМЫ................................................ Оглавление 5.2. ЛОГИЧЕСКАЯ И ФИЗИЧЕСКАЯ МОДЕЛИ ДАННЫХ........................................ 5.3. ОБОЛОЧКА...................................................................................................... 5.4. ИНТЕРФЕЙСЫ ОБМЕНА ДАННЫМИ МЕЖДУ ПРИЛОЖЕНИЯМИ................ ПРЕПРОЦЕССОР ДЛЯ ПЛОСКИХ И ОСЕСИММЕТРИЧНЫХ ЗАДАЧ ПОСТПРОЦЕССОР Г л а в а 6. ПРИМЕРЫ ЗАДАЧ, РЕШАЕМЫХ С ПОМОЩЬЮ ППП “AKORD” 6.1. ТЕСТОВЫЕ РАСЧЕТЫ.................................................................................. ТЕСТ №1. РЕШЕНИЕ УРАВНЕНИЯ ЛАПЛАСА В ЗАДАННОЙ ОБЛАСТИ ТЕСТ №2. ТЕСТИРОВАНИЕ ИТЕРАЦИОННОГО АЛГОРИТМА ТЕСТ №3. ДВИЖЕНИЕ УЕДИНЕННОЙ ВОЛНЫ ПО РОВНОМУ ДНУ ТЕСТ №4. ДВИЖЕНИЕ СОЛИТОНА НАД ПРЯМОУГОЛЬНЫМ УСТУПОМ ТЕСТ №5. ЗАДАЧА ОБТЕКАНИЯ ПРОФИЛЯ ЖУКОВСКОГО ТЕСТ №6. СХЛОПЫВАНИЕ СФЕРИЧЕСКОЙ ГАЗОВОЙ ПОЛОСТИ (ЗАДАЧА РЭЛЕЯ) ТЕСТ №7. СХЛОПЫВАНИЕ ГАЗОВОГО ПУЗЫРЯ ВОЗЛЕ ТВЕРДОЙ СТЕНКИ 6.2. ПРИМЕРЫ РЕШЕННЫХ ЗАДАЧ.................................................................... ЗАДАЧА №1. ОПРЕДЕЛЕНИЕ ФОРМЫ СВОБОДНОЙ ГРАНИЦЫ ПРИ ОБТЕКАНИИ ПОЛУКРУГОВОГО ЦИЛИНДРА ЗАДАЧА №2. ИНТЕГРАЛЬНЫЕ ХАРАКТЕРИСТИКИ УЕДИНЕННЫХ ВОЛН ЗАДАЧА №3. НАКАТ УЕДИНЕННОЙ ВОЛНЫ НА ТВЕРДУЮ СТЕНКУ ЗАДАЧА №4. ВЗАИМОДЕЙСТВИЕ СОЛИТОНОВ ЗАДАЧА №5. ДВИЖЕНИЕ СОЛИТОНА НАД ДНОМ С ПОЛУКРУГОВЫМ ВЫСТУПОМ ЗАДАЧА №6. ГОРИЗОНТАЛЬНОЕ ДВИЖЕНИЕ ПОЛУКРУГОВОГО ЦИЛИНДРА ИЗ СОСТОЯНИЯ ПОКОЯ ПО РОВНОМУ ДНУ ЗАДАЧА №7. СОВМЕСТНОЕ ВЛИЯНИЕ ТВЕРДЫХ СТЕНОК И СВОБОДНЫХ ГРАНИЦ НА ЭВОЛЮЦИЮ ГАЗОВОГО ПУЗЫРЯ ЗАДАЧА №8. ЭВОЛЮЦИЯ ПУЗЫРЯ В "СИЛЬНО" ОГРАНИЧЕННОМ ОБЪЕМЕ ЗАДАЧА №9. ИССЛЕДОВАНИЕ ЯВЛЕНИЙ НА ПОВЕРХНОСТИ ВОДЫ ПРИ СХЛОПЫВАНИИ ГАЗОВОЙ ПОЛОСТИ ЛИТЕРАТУРА ПРИЛОЖЕНИЯ ПРИЛОЖЕНИЕ 1. РАБОТА С ПРЕПРОЦЕССОРОМ ДЛЯ ПЛОСКИХ И ОСЕСИММЕТРИЧНЫХ ЗАДАЧ........................................................... ПРИЛОЖЕНИЕ 2. РАБОТА С ПОСТПРОЦЕССОРОМ....................... ПРИЛОЖЕНИЕ 3. ER-ДИАГРАММА МОДЕЛИ ДАННЫХ, ИСПОЛЬЗОВАННОЙ В ИС................................................................. ПРИЛОЖЕНИЕ 4. СТРУКТУРА ТАБЛИЦ БАЗЫ ДАННЫХ................ ПРИЛОЖЕНИЕ 5. СТРУКТУРА ФАЙЛОВ ОБМЕНА ДАННЫМИ МЕЖДУ ПРИЛОЖЕНИЯМИ ИС....................................................................... Введение ВВЕДЕНИЕ Проблемам проведения численных расчетов при решении сложных за дач математической физики с использованием ресурсов ЭВМ посвящено не мало книг и монографий. Однако немалый опыт работы авторов со студен тами дает основание утверждать, что информационных источников, где бы использовался комплексный подход при описании методов, технологий и средств численного моделирования задач на ЭВМ, по-прежнему недостаточ но. Сказываются также понятные затруднения при использовании соответст вующей литературы в учебном процессе.

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

Основное внимание уделено двум аспектам – описанию численных ме тодов граничных элементов и технологии их практического использования при решении задач с помощью пакета прикладных программ "AKORD", раз рабатываемого в течение ряда лет в Центре новых информационных техноло гий КемГУ для проведения численных расчетов в задачах гидродинамики со свободными границами.

Методы граничных элементов стали применяться сравнительно недавно.

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

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

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

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

Четвертая глава посвящена описанию структуры пакета и компонен тов, использованных при его проектировании.

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

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

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

В работе над пакетом прикладных программ активное участие принима ли аспиранты Коротков Г.Г. – препроцессор и постпроцессор для двумерных задач, Долаев Р.Р. – постпроцессор для двухмерных задач и оболочка, Стуко лов С.В. – решатель для плоских задач, выполненный на основе интегральной формулы Коши.

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

Предлагаемое учебное пособие является частью плановой работы, ко торая была определена департаментом информационных технологий Орга низации Объединенных Наций по Вопросам Образования, Науки и Культуры (UNESCO), при создании в Кемеровском государственном университете ка федры ЮНЕСКО по новым информационным технологиям в образовании и науке.

Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости ГЛАВА НЕКОТОРЫЕ СВЕДЕНИЯ ИЗ ТЕОРИИ ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ Модель идеальной несжимаемой жидкости находит широкое применение для моделирования различных гидродинамических за дач, к которым относятся и стационарные (нестационарные) задачи гидродинамики со свободными границами. Подробное описание задач гидродинамики приведено во многих монографиях и учеб никах [25, 36, 45, 47, 50, 73, 75].

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

Постановки этих задач различны, искомые функции имеют разный физический смысл. Наиболее распространенными уравне ниями этого типа являются уравнение Лапласа ( R) = 0, (1.1) или уравнение Пуассона Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости ( R) = f ( R), (1.2) где ( R ) и f ( R ) - функции радиус-вектора точки R = R ( x, y, z ) в Евклидовом пространстве.

Функцию ( R ) - зачастую можно интерпретировать как потен циал некоторого безвихревого векторного поля, описываемого уравнением V ( R) = grad ( R ). (1.3) Уравнения (1.1) и (1.2)) описывают стационарное распределе ние потенциала ( R ) внутри некоторого объема D, ограниченного поверхностью S.

Математическая задача о нахождении потенциала ( R ) форму лируется следующим образом [81].

Найти функцию ( R ), удовлетворяющую внутри D уравнению (1.2) и граничным условиям, которые могут быть заданы в одном из следующих видов:

= f1, I. на S (первая краевая задача).

= f2, II. на S (вторая краевая задача).

n + ( f3 ) = 0, III. на S (третья краевая задача);

n, f3, - это заданные функции, где f1, f - производная по n внешней нормали к поверхности S.

Первую краевую задачу для уравнения Лапласа называют «за дачей Дирихле», а вторую – «задачей Неймана». Если на одной части S1 поверхности S задано условие Дирихле, а на другой S2 =S\S1 - условие Неймана, то краевая задача называется смешан ной краевой задачей.

Если решение ищется в области D0, внутренней (или внешней) по отношению к поверхности S, то соответствующая задача назы вается внутренней (или внешней) краевой задачей.

Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости Физическое содержание каждого конкретного потенциала (при решении различных физических задач) можно найти во многих учебниках и монографиях (например, [45, 47, 50, 73]).

1.2. Основные уравнения теории идеальной жид кости Уравнение неразрывности Пусть D - фиксированная область с границей S, занятая иде альной несжимаемой и однородной жидкостью, то есть жидкостью с постоянной во всей области течения плотностью = const.

Скорость изменения массы жидкости в области D, обусловлен ная движением жидкости через границу S со скоростью v, равна n vds.

S Здесь n - вектор единичной внешней нормали к элементу поверх ности ds. Символ “x” обозначает векторное произведение, а символ " " далее по тексту означает скалярное произведение векторов.

С другой стороны, скорость изменения массы жидкости в об ласти D (элемент объема которой обозначается dV) дается форму лой с t dV.

D Поскольку материя не исчезает никуда и не возникает ниотку да, то закон сохранения массы можно выразить уравнением t dV = n vds. (1.4) D C Знак минус перед интегралом в правой части взят потому, что этот интеграл положителен, если через поверхность S за единицу вре мени вытекает больше жидкости, чем втекает, что способствует Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости уменьшению плотности во времени. Иными словами, интегралы, входящие в выражение (1.4), всегда имеют разные знаки.

Применяя к преобразованию интеграла по границе формулу Гаусса, приведем уравнение (1.4) к виду t + v dV = 0.

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

+ v = 0. (1.5) t Вводя в рассмотрение материальную производную d = + v, (1.6) dt t уравнение неразрывности (1.5) можно преобразовать к форме d + v = 0, dt которая для несжимаемой однородной жидкости с постоянной плотностью приводится к соотношению v = 0. (1.7) Если в выражении (1.3) скалярную функцию ( R ) считать по тенциалом поля скоростей, то функция, стоящая в левой части это го выражения будет представлять вектор скорости частиц жидко сти. Подставляя (1.3) в (1.7), получим уравнение Лапласа (1.1), описывающее безвихревое течение идеальной несжимаемой жид кости.

Возможность безвихревого течения идеальной жидкости дока зывает теорема Лагранжа [45].

Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости Уравнение движения Основным динамическим уравнением движения материальной точки является второй закон Ньютона [73] dv dmv F =m =, (1.8) dt dt где m – масса материальной точки, а v – ее скорость.

Произведение массы на скорость mv называется количеством движения точки.

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

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

dv dV = p nds + BdV, (1.9) dt D S D где р – давление, B – внешние массовые силы.

Если массовые силы B имеют потенциал G (т.е. B =gradG), то преобразуя в уравнении (1.9) интеграл по поверхности с помощью теоремы Гаусса в объемный интеграл, получим dv dV = ( G + p )dV, dt D D откуда следует, что dv p = grad (G + ). (1.10) dt Глава 1. Некоторые сведения из теории идеальной несжимаемой жидкости Интегралы движения Уравнение (1.10) записано в предположении, что жидкость яв ляется идеальной и несжимаемой. Если дополнительно предполо жить, что течение безвихревое, массовые силы определяются силой тяжести с вектором g, направленным противоположно оси z, то из (1.10) с учетом (1.6) следует интеграл Коши-Лагранжа:

1 p + ( ) + gz + = c(t ).

(1.11) t 2 Здесь с(t) - функция от времени, равная значению левой части в не которой точке пространства. Если жидкость на бесконечности по коится и давление на уровне z=0 равно нулю, то c(t) = 0.

Для численных расчетов нестационарных течений удобно пользоваться интегралом движения, записанным для фиксирован ной частицы жидкости. Из выражения (1.6) следует, что d + ( ) 2 и интеграл (1.11) преобразуется к виду:

= dt t d 1 p ( ) + gz + = c(t ).

(1.12) dt В случае стационарного течения, когда в фиксированной точке пространства все функции с течением времени сохраняют постоян ное значение, t = 0 и из (1.11) следует интеграл Бернулли 1 p ( )2 + gz + = c. (1.13) Глава 2. Постановка задач со свободными границами Глава ПОСТАНОВКА ЗАДАЧ СО СВОБОДНЫМИ ГРАНИЦАМИ, РЕШАЕМЫХ ПРИ ПОМОЩИ ПАКЕТА “AKORD” Основой решения всех краевых задач и задач со свободными границами, в частности, является корректность задания начальных и граничных условий или, другими словами, – корректности поста новки краевой задачи. Немаловажным является вопрос выбора ха рактерных размерных параметров для обезразмеривания постав ленной краевой задачи.

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

Стационарная задача Рассмотрим некоторую гипотетическую область решения D (рис. 1), ограниченную свободной поверхностью C1, твердой стен кой C3 и двумя границами С2 и С4, через которые протекает жид кость. Глубина жидкости на бесконечности принимается равной H, скорость течения V, давление на свободной поверхности предпо лагается равным нулю ( p = 0 ), ускорение силы тяжести g направ Глава 2. Постановка задач со свободными границами лено против направления оси y. Границы C2 и C4 должны быть вы браны на достаточно большом удалении от выступа с тем, чтобы поток через эти участки границы был близок к параллельному.

Плоский случай g Рис. 1. Схема области течения и основные обозначения Потенциал поля скоростей ( x ) внутри этой области должен удовлетворять уравнению Лапласа (1.1). На твердой границе вы полняется условие «непротекания» (условие Неймана) = 0, x C3. На боковых участках границы должны также вы n полняться краевые условия второго рода (условия втека ния/вытекания жидкости):

= V, x C2, C4. (2.1) n Свободная граница C1 является линией тока и на ней должны x C1 (t ) и дина = 0, выполняться кинематическое условие n мическое условие при р=0. Из уравнения Бернулли имеем:

1 V + gy = + gH, x C1, (2.2) 2 здесь x - радиус-вектор точки области решения, n - вектор внеш ней нормали.

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

Схематически область течения представлена на рис. 2.

g Н Рис. 2. Схема области течения и основные обозначения Пусть в расчетной области течения D(t ), ограниченной сво бодной поверхностью C1 и твердыми стенками C2 (рис. 2), решает ся нестационарная краевая задача. Если потенциал поля скоростей зависит от времени ( x, t ), то задача удовлетворяет уравнению Ла пласа. На твердой стенке C2 выполняется условие Неймана ( n = 0 ), а на свободной поверхности C1 должны выполняться кинематическое и динамическое условия:

dx =, x C1. (2.3) dt d 1 + gy = gH, x C1. (2.4) dt В качестве начальных условий необходимо задать начальное поло жение свободной границы C10 и распределение потенциала 0 на ней.

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

Осесимметричный случай Нестационарная задача Само название «осесимметричная задача» говорит о том, что у задачи существует осевая симметрия. Часто за ось симметрии при Глава 2. Постановка задач со свободными границами нимается координата z, а система декартовых координат преобра зуется к цилиндрической системе по следующему закону ( x, y, z ) (r cos, r sin, z ) и предполагается, что процесс не зави сит от угла поворота.

В новой системе координат уравнение Лапласа запишется в виде:

2 2 + + = 0. (2.5) r r 2 z r На оси симметрии r = 0 ставится условие непротекания:

= 0. (2.6) n r = Рис. 3. Переход к цилиндрической системе координат Остальные граничные условия и условия на свободной грани це остаются без изменения. Например, в случае, изображенном на (рис.3), граничные условия будут иметь вид:

d 1 + gz = 0, ( z, r ) C1 (t ), dt dr dz = =, ( z, r ) C1 (t ),, dt r dt z = 0, ( z, r ) C2.

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

Глава 2. Постановка задач со свободными границами Поскольку при помощи пакета “AKORD” решаются задачи ди намики газового пузыря в идеальной жидкости, то целесообразно привести постановку данной задачи в более общем виде [63].

Пусть в пространстве R 3 точек x =(x,y,z) рассматривается об ласть D(t), зависящая от времени t 0, ограниченная гладкими по верхностями S(t) и Г. Поверхность Г является твердой стенкой, а поверхность S(t) является границей односвязной ограниченной об ласти Q(t) (пузырь), расположенной вблизи Г (рис. 4).

z y n в b пузырь x Q(t) S(t) H Г D(t) стенка Рис. 4. Начальная геометрия пузыря возле наклонной стенки Сделаем некоторые предположения относительно движения идеальной несжимаемой жидкости. На бесконечности жидкость покоится. На жидкость действует поле сил тяжести. Давление на бесконечности постоянно и равно pa. Давление ps (t ) в пузыре Q(t) меняется по адиабатическому закону ps = C Q(t ), где Q(t ) объем области Q(t), С=const, = const 1. Массой газа заключен ного в пузыре Q(t) и силами поверхностного натяжения можно пре небречь. Начальное состояние при t=0 задается условиями: S(0) – сфера радиуса R с центром в точке (0,0,0), расстояние до поверхно сти Г равно Н, причем H R, давление в Q(0) задано и равно p0.

Краевая задача для уравнения Лапласа ставится обычным обра зом, а потенциал скорости ( x, t ) и гидродинамическое давление p( x, t ) должны удовлетворять интегралу Коши – Лагранжа d 1 1 + gz + ps = pa, x S (t ). (2.7) dt Глава 2. Постановка задач со свободными границами ps = p0 Q(t ) / Q(0) В силу сделанных предположений, Q(0) = R3.

2.2. Комплексные переменные. Плоский случай Для решения плоских краевых задач весьма эффективным яв ляется аппарат теории функций комплексного переменного, хоро шо зарекомендовавший себя в различных областях знаний [36, 45, 50].

Стационарная задача Пусть в расчетной области течения D(t), ограниченной свобод ной поверхностью C1 (t ) и твердыми стенками C2 (рис. 1), решается уравнение Лапласа:

w( z ) = 0, z = z ( x, y ) D(t ), (2.8) для функции комплексного потенциала w( z ) = ( x, y ) + i ( x, y ), где ( x, t ) - функция потенциала скоростей, ( x, y ) - функция тока, удовлетворяющие условиям Коши-Римана [45]:

= =.

;

x y y x На твердых границах выполняется условие непротекания:

= 0, z C2. (2.9) На боковых участках границы области и на дне выполняются крае вые условия вида ( x, y ) = y, z C2, C4, (2.10) ( x, y ) = 0, z C3.

Свободная граница является линией тока ( ( x, y ) = H ), на ко торой справедливо уравнение Бернулли (1.13), в данном случае имеющее вид:

Глава 2. Постановка задач со свободными границами 2 1 dw V + gy = + gH, z C1 (t ). (2.11) 2 dz Стационарная задача циркуляционного обтекания тел Задачам циркуляционного обтекания профилей посвящено множество работ, но, как правило, в них область течения безгра нична, либо ограничена твердыми стенками. Наличие свободной границы привносит дополнительную трудность в численное моде лирование данной задачи.

Пусть одиночный профиль с границей C5 обтекается потоком весомой жидкости, ограниченной свободной поверхностью C1 (t ), прямолинейным дном C3, участками втекания C2 и вытекания C (рис.5). Введем обозначения: V - скорость втекающего потока, H - глубина жидкости на бесконечности, g - ускорение силы тяжести.

Границы C2 и C4 должны быть выбраны на достаточно большом удалении от обтекаемого объекта с тем, чтобы линии тока жидко сти, втекающей через эти участки границы, были параллельны.

Для решения задачи циркуляционного обтекания профилей применять метод комплексных граничных элементов с использова нием интегральной формулы Коши для функции комплексного по тенциала w( z ) затруднительно, так как потенциал скорости ( x, y ) неоднозначно определен и при циркуляции отличной от нуля тер пит разрыв первого рода в острой кромке профиля.

Рис. 5. Схема области течения и основные обозначения Поскольку поле скоростей остается непрерывным, то данную задачу проще решать в терминах комплексно-сопряженной скоро сти W ( z) = +i = i = Vx iVy. (2.12) x x x y Глава 2. Постановка задач со свободными границами Пусть - точка, принадлежащая контуру, тогда выражение (2.12) может быть записано в виде W ( ) = Vx ( ) iVy ( ) = (Vn ( ) iVs ( ) ) ei ( ), где Vs и Vn касательная и нормальная компоненты вектора скоро сти в точке, - угол между вектором Vn и осью ОХ.

Задача об обтекании профиля может быть сведена к решению уравнения Лапласа W ( z ) = 0, z = x + iy D(t ), для аналитической в области течения D функции W ( z ). На профиле и на дне выполняется условие непротекания Vn = 0, z C5, C3. На боковых участках границы области ставятся условия втекания и вытекания жидкости Vn = V, z C2, C4. (2.13) Если обозначить касательные составляющие векторов скорости на верхней и нижней сторонах профиля через Vs+ и Vs, то условие Жуковского-Чаплыгина в острой кромке можно записать следую щим образом:

Vs+ + Vs = 0. (2.14) В силу того, что нормальные составляющие Vn+ и Vn вектора скорости на профиле равны нулю, а касательные составляющие Vs+ и Vs равны по величине и противоположны по знаку, то если вве сти обозначение Vs0 = Vs+ + Vs, условие Жуковского-Чаплыгина в острой кромке запишется следующим образом:

Vs0 = 0. (2.15) Глава 2. Постановка задач со свободными границами Свободная граница является линией тока, на которой справед ливо уравнение Бернулли V 1 W ( z ) + gy = + gH, z C1. (2.16) 2 Нестационарная задача Пусть в расчетной области D(t) (рис. 5) решается нестационар ная задача для уравнения Лапласа (2.8). Краевые условия на твер дой границе имеют вид (2.9), а на свободной поверхности выпол няются кинематическое и динамическое условия:

dz = +i z C1 (t ). (2.17), dt x y d 1 dw + g ( y H ) = 0, z C1 (t ). (2.18) dt 2 dz Кроме того, необходимо задать начальные условия - положение свободной границы C10 и распределение потенциала 0 на ней.

2.3. Приведение краевых задач к безразмерному виду При проведении численных расчетов целесообразно перейти в исходных уравнениях к безразмерным переменным. Для этого не обходимо выбрать некоторые характерные для решаемой задачи независимые размерные переменные. Подробную информацию о методах обезразмеривания можно найти во многих учебниках по гидродинамике, например [45, 73].

Стационарные задачи Для стационарных задач наиболее часто выбираются в качестве независимых размерных величин скорость втекающего потока V и глубина канала H.

x = x / H, y = y/H, Введем безразмерные переменные = / HV, V = V/V и сделаем замену переменных в уравнениях, описывающих стационарную краевую задачу.

Глава 2. Постановка задач со свободными границами Все уравнения и граничные условия останутся без изменения за исключением условий (2.1), (2.2), (2.11) и (2.13), которые преоб разуется к виду В вещественных переменных В комплексных переменных Vn = 1, Vn = 1, 1 y 1 y dw = 1 + 2, = 1+ 2 2, (2.19) Fr 2 dz Fr V где Fr = - число Фруда, (черта над безразмерными перемен gH ными опущена).

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

Безразмерные переменные в этом случае могут быть введены следующим образом:

x = x / H, y = y/H, = / H gH, V = V/ gH, t = t/ H / g, p = p / gH.

Краевая задача и граничные условия останутся без изменений, за исключением уравнений (2.4) и (2.18).

В вещественных переменных В комплексных переменных d 1 d 1 dw + ( y H ) = 0, + ( y H ) = 0, (2.20) dt 2 dt 2 dz здесь, как и раньше черта над безразмерными переменными опу щена.

Глава 2. Постановка задач со свободными границами Когда в нестационарных задачах речь идет об эволюции газо вого пузыря, то характерными размерными величинами, как прави ло, являются максимальный радиус пузыря R и ускорение силы тя жести g. Безразмерные переменные записываются аналогично пре дыдущему случаю.

Как и раньше все исходные уравнения остаются без изменений, и лишь интеграл Коши-Лагранжа (2.7) преобразуется к виду:

d 1 + z = + q (t ), dt где p pa ps =, = s, q (t ) = 1 Q (t ) / Q(0).

Rg Rg Глава 3. Математические методы и алгоритмы Глава МАТЕМАТИЧЕСКИЕ МЕТОДЫ И АЛГОРИТМЫ, ИСПОЛЬЗУЕМЫЕ В ППП «AKORD»

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

Свести задачу для дифференциального уравнения эллиптиче ского типа к граничному интегральному уравнению можно различ ными способами. Ниже описаны методы, реализованные в настоя щей версии пакета прикладных программ «AKORD» - комплексный метод граничных элементов (КМГЭ), построенный на использова нии интегральной формулы Коши [32], и метод граничных элемен тов (МГЭ), построенный на использовании третьей формулы Грина [21, 24, 78].

3.1. КМГЭ на основе интегральной формулы Коши Для области D, ограниченной кусочно-гладкой замкнутой гра ницей C = C1 C2, справедлива интегральная формула Коши, кото Глава 3. Математические методы и алгоритмы рую можно записать с помощью предельных формул Сохоцкого в виде:

w( z ) w( z0 ) = dz, (3.1) ( z0 )i C z z где ( z0 ) = 2 для внутренней точки, ( z0 ) = для точки на глад кой границе C, ( z0 ) = для угловой точки границы C ( - угол при вершине). Положительное направление обхода контура C бе рется таким образом, чтобы область D оставалась слева.

Учитывая, что на свободной поверхности известна действи тельная часть ( x, y ) функции w( z ), а мнимая часть ( x, y ) из вестна на твердых стенках, мы имеем смешанную краевую задачу для функции w( z ).

Численное решение задачи можно получить, разбив контур C на N линейных элементов j с узлами z j ( j = 1, N ). Тогда w( z ) = lim G ( z ), где G(z) - линейная глобальная пробная функ max j N N ция для z j и G ( z ) = w j j ( z ), w j - значение w( z ) в точке j =1 j = z j, j ( z ) - линейная базисная функция:

( z z j ) /( z j z j 1 ), z j 1, j ( z ) = ( z j +1 z ) /( z j +1 z j ), z j, 0, z j 1 j.

Вычисление интегралов После указанного разбиения и линейной аппроксимации функ ции w( z ) на границе интеграл Коши можно вычислить аналитиче ски в смысле главного значения при z z j.

G ( ) 2 iw( z j ) = lim d.

zz j z C Глава 3. Математические методы и алгоритмы Обозначение z z означает, что точка z стремится к точке z j, j оставаясь все время внутри области D.

G ( ) G ( ) lim d = lim d + zz j z z zz j C j. (3.2) G ( ) G ( ) N d + d + lim z zj zz j m = j +1 m m j m j + G ( ) d. В выражении (3.2) интегра Введем обозначение: I m = zj m лы под знаком суммы не имеют особенности и вычисляются сле дующим образом:

zm + G ( ) wm m + wm +1 m + d = d = Im = zj zj m zm (3.3) ( z j zm ) wm +1 ( z j zm +1 ) wm zm +1 z j wm +1 wm + ln.

zm +1 zm zm z j zm +1 zm Два оставшихся интеграла имеют особенности при z z j и могут быть представлены в виде:

z j + zj w j 1 j 1 + w j j w j j + w j +1 j + d + lim d = I j + I j +1 = lim zj zj zz j zz j z j 1 zj w w j d z zj z z j 1 z zj z j = lim z d + zi zi 1 w j zi zi 1 w j 1 z + j z z j z j z j 1 z j z j Глава 3. Математические методы и алгоритмы w w z j +1 z j + z zj z z j + z d j + + lim j d + w j +1 wj.

zi +1 z j z z z z j z j +1 z j zi +1 zi zj zj Комбинируя четные и нечетные слагаемые по индексу j, по лучим:

z j +1 z j I j + I j +1 = w j +1 w j 1 + w j ln z j 1 z j.

Комплексный логарифм в последнем слагаемом расписывается на действительную и мнимую части следующим образом:

z j +1 z j z j +1 z j + i arg( z j +1 z j ) arg( z j 1 z j ) = = ln ln z j 1 z j z j 1 z j z j +1 z j + i (2 ( z j )), = ln z j 1 z j где ( z j ) - угол при вершине, образованный j и j+1 элементами.

Таким образом, выражение (3.2) окончательно может быть за писано:

z j +1 z j N iw j = w j +1 w j 1 + w j ln + Im, (3.4) z j 1 z j m = m j, j + где интегралы I m вычислены по формуле (3.3).

Подставив в равенство (3.4) известные действительные или мнимые части функции w( z ) при j=1,2,…,N, получим систему N N линейных алгебраических уравнений BX=F, (3.5) Глава 3. Математические методы и алгоритмы здесь В - полнозаполненная матрица, X - вектор неизвестных для определения Re w на C2 и Im w на C1, F - вектор правой части.

Для решения системы (3.5) используется метод Гаусса с выбо ром ведущего элемента [69, 84].

3.2. Модифицированный КМГЭ В ходе решения задачи о циркуляционном обтекании тел на границах области известны касательные и нормальные компоненты вектора скорости Vs и Vn. Поэтому интегральную формулу Коши с помощью предельных формул Сохоцкого для аналитической в об ласти D функции комплексно-сопряженной скорости W ( z ) можно записать в следующем виде:

W ( z0 ) = (Vn ( z0 ) iVs ( z0 ) ) ei ( z0 ) = (Vn ( z ) iVs ( z ) ) ei ( z ) dz, (3.6) = ( z0 )i C z z где ( z0 ) = 2 для внутренней точки, ( z0 ) = для точки на глад кой границе C, ( z0 ) = для угловой точки границы C ( - угол при вершине). Положительное направление обхода контура C бе рется таким образом, чтобы область D оставалась слева, (внешняя граница обходится против часовой стрелки, внутренняя - по часо вой). Зная Vs и Vn, можно найти Vx и V y по формулам:

Vx = Vn cos Vs sin, V y = Vs cos + Vn sin. (3.7) Поскольку в ходе итерационного процесса на свободной гра нице известна касательная составляющая скорости, а на дне, боко вых стенках и профиле - нормальная составляющая, то для функ ции W ( z ) имеем смешанную краевую задачу.

Разобьем границу области на N линейных элементов j с уз лами z j. Воспользовавшись уравнением (3.4), записанным для функции комплексно-сопряженной скорости W ( z ), расписанной по Глава 3. Математические методы и алгоритмы касательным Vs и нормальным Vn составляющим вектора скорости, получим:

( ) ( ) ( ) i j i j +1 i j i Vn j iVs j e = Vn j +1 iVs j +1 e Vn j 1 iVs j 1 e + z j +1 z j ( ) N i j + Vn j iVs j e + Im, (3.8) ln z j 1 z j m = m j, j + где ( ) ( ) I m = Vnm+1 iVsm+1 e i m+1 Vnm iVsm e i m + ( ) ( ) e i m+1 ( z j zm +1 ) Vnm iVsm e i m ( z j zm ) Vn iVs z z ln m +1 j + m +1 m + zm z j zm +1 zm zm +1 zm, z z z j z j 1 z j +1 z j z j z j i j = i j +1 j + +.

e z z 2 z z 2 2 z j +1 z j z j z j j +1 j j j Записывая уравнение (3.8) для каждого узла границы и разде ляя мнимые и действительные части, получим AX+iBX=0, где A и B - полнозаполненные матрицы N 2 N ( N - строк, 2N - столб цов), вектор X=X(Vs1,Vn1,Vs2,Vn2,...,VsN,VnN ). Далее, следуя [32], по лучаем следующую систему уравнений GY=F, (3.9) где матрица G и вектор правой части получаются следующим обра зом: если в узле z j задана касательная составляющая скорости Vs j, то берется j строка матрицы B и после выборки элементов строки, соответствующих неизвестным значениям Vs или Vn во всех ос тальных узлах, получается j строка матрицы С, j элемент вектора Y будет Vn j, j элемент вектора F - сумма известных значений Vs Глава 3. Математические методы и алгоритмы или Vn, умноженных на соответствующие элементы матрицы B;

ес ли в узле z j задана нормальная составляющая скорости Vn j, то для построения j строки системы (3.9) используется матрица A.

Случай безграничной области Интегральная формула Коши справедлива и для безграничной области D, если функция W ( z ) обращается в нуль на бесконечно сти. При обтекании профиля безграничным потоком жидкости комплексный потенциал на бесконечности имеет следующее раз ложение:

w( z ) = ( x, y ) + i ( x, y ) = K + V ei z + ln z, 2 i где K - произвольная комплексная постоянная, Г - циркуляция ско рости вдоль контура профиля, - угол атаки профиля, V - модуль вектора скорости потока на бесконечности.

Комплексно-сопряженная скорость потока на бесконечности W ( z ) = Vx ( x, y ) iVy ( x, y ) = V ei. (3.10) Уравнение (3.6) с учетом (3.10) примет вид:

(Vn ( z ) iVs ( z ) ) ei ( z ) dz.

i W ( z0 ) = V e + ( z0 )i C z z В данном случае контур обходится по часовой стрелке. При по строении системы уравнений (3.9) следует добавить слагаемое в известный правый вектор F (мнимую или действительную часть выражения V ei к каждому элементу вектора в зависимости от заданного граничного условия в j узле).

3.3. МГЭ на основе третьей формулы Грина Плоский случай Поскольку потенциал течения идеальной несжимаемой жидко сти является гармонической функцией, можно воспользоваться из Глава 3. Математические методы и алгоритмы вестной из теории гармонических функций для уравнения Лапласа функцией Грина, которая для двумерных задач имеет вид:

G (, x) = ln(r ), где = ( xi, yi ) - точка наблюдателя, x=(x,y) – текущая точка грани цы, - расстояние между точками и x:

r r (, x) = [( xi x) 2 + ( yi y ) 2 ]1/ 2.

Используя третью формулу Грина, можно записать следующее интегральное уравнение:

G (, x) ( x) ( ) ( ) + ( x) d ( x) = G (, x)d ( x), (3.11) n( x ) n( x ) где ( x) = 1 2 - граница области D, ( x) - гармоническая функ ция, n ( x) - внешняя по отношению к области D единичная нормаль к поверхности ( x). Параметр ( ) определяется следующим об разом: ( ) = 2 для внутренней точки, ( ) = для точки на глад кой границе, ( ) = для угловой точки границы ( - угол при вершине).

Как видно из (3.11) гармоническая функция определяется в произвольной точке области по значениям самой функции и ее нормальной производной n, известным на границе.

Обычно для краевых задач на границе Г задается либо функция (задача Дирихле), либо функция / n (задача Неймана), либо на одной части границы ( Г1 ) - функция, а на другой ( Г 2 ) - функ ция / n (смешанная краевая задача). Во всех случаях на гра нице области имеет место граничное интегральное уравнение (3.11), решение которого в аналитическом виде может быть полу чено только для очень узкого круга задач. Поэтому для нахожде ния решения применяются численные методы.

Численное решение интегрального соотношения (3.11) можно получить, разбив границу Г на N линейных элементов Г j узлами Глава 3. Математические методы и алгоритмы ( x j, y j )( j = 1, N ) так, чтобы j - му элементу принадлежали узлы ( x j, y j ) и ( x j +1, y j +1 ).

Введем глобальную базисную функцию N 2j 1 ( x, y ), ( x, y ) Г j 1, N1j ( x, y ), ( x, y ) Г j, j ( x, y ) = 1, ( x, y ) = ( x j, y j ), 0, ( x, y ) Г j 1 Г j, где N k,(k = 1, 2) - локальные интерполирующие функции на эле j Г j от однородной координаты менте такой, что = (2 s / L j 1) [ 1,1], ( L j - длина j - го элемента, s (0, L j ) - те кущая координата). Координата s может быть представлена в виде:

s = ( x x j +1 ) 2 + ( y y j +1 ) 2, ( x, y ) Г j.

На j - м элементе локальные интерполирующие функции пред ставляются как:

1 1+ N1j = N 2j =. (3.12) ;

2 Введем, как и в п.3.1, пробную функцию N R ( x, y ) = f j j ( x, y ), (3.13) j = где f j = f ( x j, y j ) - значение непрерывной функции f ( x, y ) в j -ом узле.

Очевидно, что f ( x, y ) = lim R ( x, y ).

max j Глава 3. Математические методы и алгоритмы Интегральное уравнение (3.11) после подстановки в него функ ций (3.13) может быть приведено к следующему виду:

N N i + H ij j + H iii = Gij q j + Gii qi, (3.14) j =1 j = i j i j где несингулярные интегралы имеют вид:

j qi (s)ds;

ji (s)ds, * * H ij = Gij = Pj Pj здесь Pj = j 1 j - носитель глобальной базисной функции j.

Сингулярные интегралы 3 H ii 0;

Gii = Li 1 ln Li 1 + Li ln Li, 4 2 2 где L j - длина j-го элемента.

Приведение интегралов Входящие в (3.14) интегралы могут быть приведены к виду, удобному для дальнейших вычислений. Имеем:

N 2j 1qij 1 ( s )ds + j qi ds = N1 qij (s)ds = * * j* H ij = j 1 j Pj 1 (1 + ) (1 ) D j 1 Dj r 2 d = I1 + I 2, d = 16 1 rij21 16 ij где 1 1 1 d d d d D j 1 Dj D j 1 Dj r2 ;

r2 I1 = I2 = +, 16 1 rij21 16 16 16 1 rij 1 ij 1 ij ( )( y j +1 y j ) ( y j + y j +1 2 yi )( x j +1 x j ), D j = x j + x j +1 2 xi Глава 3. Математические методы и алгоритмы ( ) + ( y j + y j +1 2 yi + ( y j +1 y j ) ) rij2 = x + x j +1 2 xi + ( x j +1 x j ) j.

Аналогично преобразуются интегралы Gij.

N 2j 1ij 1 ( s )ds + jij (s)ds = N1 ij (s)ds = * * j* Gij = j 1 j Pj 1 L j 1 Lj (1 )ln rij d = I3 + I 4, (1 + )ln rij21d = 16 1 где 1 L j 1 Lj ln rij d ;

ln rij21d I3 = 16 1 1 L j 1 Lj ln rij d.

ln rij21d I4 = + 16 1 Аналитическое вычисление интегралов Интегралы I k,(k = 1, 4) состоят из двух слагаемых I k = I kj 1 + I kj и отличаются друг от друга только индексами. Поэтому проведем дальнейшие вычисления интегралов только для j -го индекса:

1 1 1 d d r2 ;

ln r = ln r 2 d.

I1j I 2j I 3j I 4j = = = ;

d ;

r 1 1 1 Введем дополнительные обозначения:

r 2 = ( A 2 + 2 B + C );

A = L2 ;

L = Lj;

B = ( x j +1 x j )( x j +1 + x j 2 xi ) + ( y j +1 y j )( y j +1 + y j 2 yi );

Глава 3. Математические методы и алгоритмы C = ( x j +1 + x j 2 xi ) 2 + ( y j +1 + y j 2 yi ) 2 ;

D = D j =B2 -AC.

Геометрическая интерпретация величин A, B, C, D Рассмотрим элемент j границы области C и расположим точку наблюдения ( xi, yi ) как показано на рис.6.

Рис. 6. Пояснительная схема к аналитическому вычислению интегралов Введем вектора x j +1 + x j y j +1 + y j {( x )}.

)( P = xi, yi, F = x j, y j +1 y j j + 2 Тогда коэффициенты A, B, C можно записать через скалярные про изведения A = ( F F );

B = 2( P F );

C = ( P P), отсюда 2 B 2 AC = 4 P F (cos 2 ( P, F ) 1) 0. (3.15) Выражение (3.15) обращается в нуль, если вектора P и F кол линеарны. Это может быть только в случае, когда точка наблюда теля либо лежит на рассматриваемом элементе и совпадает с одним из узлов (в этом случае D 0, A C = ± B, знак “+” или ”-” зависит от того, с каким узлом совпадает точка наблюдателя), либо точка наблюдателя лежит на прямой, проходящей через элемент j.

Вычисление интегралов для случая B 2 AC Глава 3. Математические методы и алгоритмы Проводя несложные выкладки, можем получить аналитические выражения для интегралов I k,(k = 1, 4) B+ A B A I1 = arctg arctg, CA B 2 CA B 2 CA B 2 2M B I 2 = ln 1 I1, A M2 A B C I 3 = 4(ln 2 + 1) + ln( M1 M 2 ) + I 2 + I1, 2 1 C M1 2 B AC B I 4 = 1 ln + + I2, 2 A M2 A 2A M1 = A + 2 B + C ;

M 2 = A 2B + C.

где Вычисление интегралов для случая B 2 AC = 0, но A C В этом случае интегралы примут вид:

A 2A 2B I1 =, I2 = 2 + 2 arth, B 2 A2 B A2 A B 1 M 2 B 2C A A 4B arth, I 4 = ln 1 + I 3 = 4 + ln( M1M 2 ) + arth.

A B 2 M2 A A B Численное интегрирование Интегралы I k,(k = 1, 4) могут быть найдены численно с помо щью стандартных квадратурных формул Гаусса [46].

1 d d n n r 2 = Ai / r (i ), = Aii / r 2 (i ), r i =1 i = 1 1 n n d = Ai ln r (i ), d = Aii ln r 2 (i ), ln r ln r 2 2 i =1 i = 1 Глава 3. Математические методы и алгоритмы где Ai,i - веса и абсциссы квадратурной формулы Гаусса.

Квадратурные формулы Гаусса Формулы Гаусса являются наиболее распространенными фор мулами численного интегрирования, как формулы, имеющие наи высший порядок точности при заданном числе узлов в области ин тегрирования.

Таблица Веса и абсциссы квадратурной формулы Гаусса Абсциссы i N Веса Аi 1.00000000000000 0. 1.00000000000000 -0. 0.888888888888889 0. 0.555555555555556 0. 3 0.555555555555556 -0. 0.652145154862546 0. 0.652145154862546 -0. 0.347854845137454 0. 4 0.347854845137454 -0. 0.568888888888889 0. 0.478628670499366 0. 0.478628670499366 -0. 5 0.236926885056189 0. 0.236926885056189 -0. 0.417959183673469 0. 0.381830050505119 0. 0.381830050505119 -0. 0.279705391489277 0. 7 0.279705391489277 -0. 0.129484966168870 0. 0.129484966168870 -0. В одномерном случае погрешность квадратуры Гаусса для вычисления интеграла 1 n f ( )d Ai f (i ) +, J= i = имеет оценку 22 n +1 (n!) 4 d 2n f.

max 2n + 1 ((2n)!)3 [ 1,1] d 2 n Глава 3. Математические методы и алгоритмы Квадратурные формулы Гаусса точны для многочленов степе ни 2n-1. Значения весов Ai и абсцисс i приведены в таблице 1.

Осесимметричный случай При осесимметричном потенциальном течении в цилиндриче ских координатах (r, z, ) функции и n не зависят от угло вой координаты. Поэтому интегральное уравнение (3.11) может быть преобразовано к виду:

G(, x) ( ) ( ) + ( x) r ( x)ds ( x) = n Г (3.16) ( x) = G(, x)r ( x)ds ( x).

n Г Фундаментальное решение в цилиндрических координатах имеет вид:

d G(, x) =, 4 a b cos где a = r 2 ( ) + r 2 ( x) + [ z ( ) z ( x) ] 2 ;

b = 2r ( )r ( x) ;

r ( x), z ( x) и r ( ), z ( ) - координаты фиксированной и переменной точек в ме ридиональной плоскости;

- контур пересечения тела вращения с меридиональной плоскостью, n - внешняя нормаль по отношению к области интегрирования.

Интеграл для функции G(, x) сводится к полному эллиптиче скому интегралу первого рода. Для этого следует интервал интег рирования разбить на два подинтервала (0, ) и (,2 ) и перейти к новым переменным = 2 и = + 2. В результате получа ем:

G(, x) = K (m) / (a + b)1/ 2, (3.17) Глава 3. Математические методы и алгоритмы n [ E ( m) K ( m) ] G(, x) G G nz = r = nr + n r z 2 r ( x) a + b nr [ r ( ) r ( x) ] + nz [ z ( ) z ( x) ] E (m), ( a b) a + b 4r ( )r ( x) m2 = 0 m 1.

, [ r ( ) + r ( x)] + [ z ( ) z ( x)] 2 Здесь nr, nz - компоненты нормального вектора в направлении осей r и z, соответственно;

K (m) и E (m) - полные эллиптические интегралы, для вычисления которых можно использовать их поли номиальные представления из работы [1]:

n E (m) = 1 + [ci m1 di m1 ln(m1 )] + (m);

i i i = n K (m) = [ai m1 bi m1 ln(m1 )] + (m);

i i (3.18) i = [r ( x) r ( )]2 + [ z ( x) z ( )] m1 = 1 m =.

2 [r ( x) + r ( )] + [ z ( x) z ( )] При n = 4 остаточный член удовлетворяет оценке (m) 2 108, а ai, bi, ci и di - константы, представленные в таблице 2.

Необходимо заметить, что в осесимметричном случае функция Грина зависит не только от расстояния между двумя точками и r, но и от расстояния между рассматриваемой точкой и осью симметрии.

Таблица Коэффициенты в разложении эллиптических интегралов a0 =1.38629436112 b0 =0.5 c0 =0.0 d0 =0. a1 =0.09666344259 b1 =0.12498593597 c1 =0.44325141463 d1 =0. a2 =0.03590092383 b2 =0.06880248576 c2 =0.06260601220 d 2 =0. a3 =0.03742563713 b3 =0.03328355346 c3 =0.04757383546 d3 =0. Глава 3. Математические методы и алгоритмы a4 =0.01451196212 b4 =0.00441787012 c4 =0.01736506451 d 4 =0. Таблица Веса и абсциссы пятиточечной квадратуры Гаусса с логарифмическим весом Абсциссы i Веса Аi 0.29789 34717 82894 0.02913 44721 0.34977 62265 13242 0.17397 72133 0.23448 82900 44052 0.41170 25202 0.09893 04595 16633 0.67731 41745 0.01891 15521 43196 0.89477 13610 Вычисление сингулярных интегралов Из выражений (3.18) видно, что все диагональные элементы матриц нii и Gii будут иметь логарифмическую особенность. Для того, чтобы применить формулы численного интегрирования к по строению диагональных элементов, необходимо сингулярный ин теграл разбить на две части - сингулярную и несингулярную. Для сингулярной части можно применить специальную квадратурную формулу Гаусса с логарифмическими весами:


f (2 n ) ( ) n J = f ( )ln(1/ )d Ai f (i ) + K. (3.19) (2n)!

i = Формула (3.19) точна, если функция f ( ) - многочлен степени не выше (2n-1). Веса и абсциссы пятиточечной квадратуры (3.19) приведены в таблице 3. Для этого случая К=5 106.

3.4. Методы решения систем линейных алгебраи ческих уравнений Решение краевых задач численными методами сопряжено с решением больших систем линейных алгебраических уравнений.

Все методы решения линейных систем можно разделить на две группы: прямые (Гаусса, Краута, Холецкого и т.д.) и итерационные (Гаусса-Зейделя, минимальных невязок, нелинейные регуляризую щие и т.д.). Выбор метода для решения конкретной системы зави сит от ее свойств. В методе граничных элементов матрицы систем уравнений получаются полностью заполненными. В случае задачи Неймана исходный эллиптический дифференциальный оператор Глава 3. Математические методы и алгоритмы сводится к интегральному уравнению Фредгольма II рода. Полу ченная для него в результате применения метода граничных эле ментов матрица системы линейных алгебраических уравнений об ладает “хорошими” свойствами за счет диагонального преоблада ния и может быть решена прямым или итерационным методом.

Если требуется решить задачу Дирихле, то в результате получается граничное интегральное уравнение Фредгольма I рода, которое яв ляется некорректным по Адамару, что, вообще говоря, приводит к “плохо” обусловленной системе линейных алгебраических уравне ний [80]. Поэтому для численного решения некорректных задач требуется применение методов, основанных на различных способах регуляризации [80]. В пакете “AKORD” используется нелинейный регуляризирующий алгоритм, предложенный в работе [82].

Эффект саморегуляризации С другой стороны, если линейное интегральное уравнение Фредгольма I рода имеет логарифмическую особенность в ядре, то есть возможность решить его напрямую сведением к системе ли нейных алгебраических уравнений, не применяя какого-либо регу ляризирующего алгоритма. Квадратурные формулы Гаусса и куба турные формулы Радо (пространственный случай) для решения уравнений Фредгольма I рода с ядрами типа функций Грина поро ждают саморегуляризирующие алгоритмы, в которых параметром регуляризации является шаг квадратурной формулы [27, 28].

Число обусловленности Важным критерием близости матрицы A к вырожденной явля ется число обусловленности cond ( A) = A A1, которое играет роль множителя в увеличении относительной ошибки. Для опреде ления числа обусловленности применяется алгоритм, изложенный в монографии [84]. Из этой же работы в ППП «AKORD» включены программы DECOMP и SOLVE, осуществляющие декомпозицию матрицы А и ее решение прямым методом Гаусса с выбором веду щего элемента [69].

Метод Гаусса – Зейделя Среди итерационных методов, применяемых в МГЭ для реше ния систем линейных алгебраических уравнений, популярным яв ляется метод Гаусса – Зейделя с релаксацией.

Пусть задана система линейных алгебраических уравнений Глава 3. Математические методы и алгоритмы n aik xk = bi, i = 1, 2,..., n.

k = Рекуррентное соотношение вида 1 i n = aik xk + aik xkj bi, xi* * aii k =1 k =i xij +1 = xij + ( xi* xij ) определяет метод Гаусса – Зейделя с релаксацией, где - параметр релаксации (0 1). Метод релаксации является обобщением итерационного метода Гаусса-Зейделя и при = 1 точно с ним сов падает.

Вычислительные затраты При реализации численных методов возникает естественный вопрос о точности метода и о вычислительных затратах. Решая плоские и осесимметричные задачи с помощью КМГЭ и МГЭ, по лучают полнозаполненные матрицы систем линейных алгебраиче ских уравнений. Число этих уравнений, как правило, менее 1000, а максимальное число обусловленности в самых худших случаях бы ло порядка 103. Решение таких систем осуществлялось методом Гаусса с выбором ведущего элемента. На решение системы в дан ном случае требуется N 3 операций.

Таблица Параметры релаксации, на пространственных сетках, близкие к оптимальным Способ разбиения Количество Параметр релаксации в Параметр релаксации при узлов на сфе- безграничной жидкости наличии твердой стенки ре /число итераций /число итераций На 4 зоны 402 0.57 / 28 0.50 / 578 0.55 / 30 0.48 / На 6 зон 386 0.59 / 28 0.52 / 866 0.54 / 33 0.47 / Икосаэдр 162 0.35 / 50 0.27 / 642 0.21 / 90 0.16 / В случае решения пространственных задач количество уравне ний системы значительно увеличивается. На одну итерацию в ите рационных методах требуется n 2 операций и, когда итерационный Глава 3. Математические методы и алгоритмы процесс сходится меньше чем за n шагов, привлекательность ите рационных методов перед прямыми налицо. При решении СЛАУ методом Гаусса – Зейделя с релаксацией необходимо задать опти мальное значение параметра opt, при котором удается получить наилучшую скорость сходимости. К сожалению, аналитическая оценка этого параметра имеется для очень узкого класса задач.

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

Это методы, основанные на зонном разбиении поверхностей и на «дроблении» икосаэдра. Значения параметра релаксации приведе ны для задачи обтекания сферы безграничным потоком жидкости и для задачи обтекания сферы вблизи твердой стенки.

3.5. Алгоритмы построения свободных границ Стационарная задача Итерационная схема Пусть на к-ом итерационном шаге известна форма свободной границы C1 (рис. 1). Для того, чтобы продолжить итерационный процесс, необходимо определить значение потенциала на грани це C1. Поскольку скорость на границе C1 направлена по касатель ной к контуру, то s = q или ( dW dz = q ).

Так как потенциал скорости определяется с точностью до ад дитивной константы, полагаем 1 = 0. Далее для любой точки сво бодной границы имеем qi + qi + i +1 = i + si, (3.20) _ где i = 1, N g 1 - номера узлов точек свободной границы, qi = q ( yi ) ( xi +1 xi )2 + ( yi +1 yi ) определяется формулой (2.19), где si = Глава 3. Математические методы и алгоритмы длина i -го элемента свободной границы, N g - количество узлов на свободной границе C1.

Определение формы свободной границы Алгоритм нахождения свободной границы осуществляется по следующей схеме:

k 1) задается некоторое положение свободной границы C1 ;

2) по формуле (3.20) определяются значения потенциала i в узлах C1k ;

3) решается смешанная краевая задача с заданными значениями по тенциала ik на свободной границе C1k ;

k 4) в точках свободной границы C1 определяются значения компо нент вектора скорости ui и vi ;

5) из условия коллинеарности вектора скорости и касательной к границе ( dy / dx = v / u ) вычисляется новое положение свободной границы C1k + yik++1 = yik +1 + yik.

Цикл повторяется до выполнения требуемой точности max yik +1 yik. (3.21) i Нахождение приращения y Существенным моментом построения свободной границы яв ляется то, что абсциссы х точек свободной границы C1 не меняют своих значений в течение всего итерационного процесса, изменяет ся только ордината yi. Будем аппроксимировать свободную грани цу C1 (на к-ой итерации) кусочным многочленом второй степени y = Ax 2 + Bx + C.

Неизвестные коэффициенты А и В на элементе могут быть найдены по значению производной в двух узловых точках свобод ной границы Глава 3. Математические методы и алгоритмы dy = 2 Ax + B. (3.22) dx Записывая правую часть этого уравнения с учетом кинематиче ского условия в узлах границы C1 и разрешая полученную систему относительно А и В, приходим к дифференциальному уравнению:

dy i i +1 x x x i i +1 i +1 i =. (3.23) dx ui ui +1 xi xi +1 ui +1 ui xi xi + Интегрируя последнее равенство в интервале [ xi, xi +1 ], получа ем новое значение ординаты на (к+1)-й итерации:

i +1 i xi +1 xi yi(+1 1) = yi( k +1) + yi( k ), k+ yi( k ) = ( + где.

) ui +1 ui Другим способом вычисления приращения yik, дающим более точные результаты, служит способ разложения в ряд Тейлора:

1 d 3 vi 1 d vi v = i ( xi +1 xi ) + ( xi +1 xi ) + ( xi +1 xi ) 2 yik + 4! dx3 ui 2! dx ui ui.

Данный итерационный алгоритм был развит в работах [18, 78].

Если речь идет о построении решений на основе функции тока или комплексного потенциала, то итерационная схема поиска сво бодной границы выглядит несколько иначе:

yik +1 = yik + yik, где yik = 1 ik, а ik - это текущее значение функции тока в i – точке свободной границы.

Приведенная к безразмерному виду краевая задача зависит от безразмерного параметра – числа Фруда Fr. Впервые о не единст Глава 3. Математические методы и алгоритмы венности возможных форм течений тяжелой жидкости при числах Фруда, близких к единице, было доказано в работе [61]. В работе [66] с помощью вариационного принципа доказано, что для беско нечного множества значений чисел Фруда задача имеет, по крайней мере, два различных решения. Следовательно, для получения одно значного решения использовать число Фруда в качестве параметра задачи нецелесообразно. В работе [35] отмечается, что данная зада ча имеет единственное решение, если в качестве определяющего течение параметра вместо числа Фруда задавать безразмерный па раметр V = V0 / V, характеризующий отношение скорости V0 в вершине волны к скорости набегающего потока V. При этом чис ло Фруда есть функция от V: Fr=Fr(V). Данная зависимость легко получается из уравнений (2.19):

y0 Fr 2 = 2. (3.24) 1V Следовательно, выражение для модуля скорости на свободной границе преобразуется к виду:

) yy 11.

( q = 1 1V 2 (3.25) Подстановка в итерационные алгоритмы формулы (3.24) вме сто уравнения (3.20) дает возможность получать два решения, на хождение которых осуществлялось различными способами: мето дом граничных элементов с использованием как функции тока, так и потенциала скоростей;


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

После того, как итерационный процесс устанавливается, число Фруда вычисляется по формуле (3.24).

Очевидно, что параметр V изменяется в пределах [0;

1): 1 - со ответствует бесконечному числу Фруда;

0 - числу Фруда, близкому к единице, при котором амплитуда построенной волны максималь на.

Глава 3. Математические методы и алгоритмы Алгоритм определения свободной границы в методе, исполь зующем комплексно-сопряженную скорость аналогичен алгоритму, изложенному выше с той лишь разницей, что в точках свободной границы ищется не потенциал i, а касательная составляющая век тора скорости vs.

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

В работе М.А. Лаврентьева и Б.В. Шабата [47] предложен ме тод решения нестационарных задач со свободными границами. Его идея состоит в редукции общей нелинейной задачи, определенной на непрерывном интервале времени, к последовательности линей ных краевых задач теории потенциала на дискретном множестве временных шагов из заданного временного интервала. Эта методи ка хорошо зарекомендовала себя и применялась в работах [4, 8, 17, 43, 91, 92].

Описанные нестационарные краевые задачи отличаются от тра диционных задач математической физики тем, что время явно входит только в граничные условия (например: (2.17) и (2.20)), представляющие собой обыкновенные дифференциальные уравне ния первого порядка, для интегрирования которых используется явный метод Эйлера.

Пусть в некоторый момент времени t0 задано положение сво бодной границы C10 и распределение потенциала 0 на ней. Далее необходимо решить уравнение Лапласа с условием Дирихле на C и условием Неймана на C2. Новое положение свободной границы и распределение потенциала на ней для момента времени t0 + мож но вычислить, используя условия (например, (2.17) и (2.20)), дис кретный аналог которых расписывается по схеме Эйлера следую щим образом:

x k +1 = x k + ( ), k Глава 3. Математические методы и алгоритмы k +1 = k + 0,5 ( ) k ( y k H ), k = 0,1, 2,..., (3.26) где х k = ( x k, y k ), k - значения функций на k-м шаге по времени.

Таким образом, для момента времени t0 + опять получаем смешанную краевую задачу для уравнения Лапласа с граничными условиями Дирихле 3.26 и условиями Неймана, соответствующими решаемой задаче. Повторное ее решение и использование гранич ных условий позволяет определить положение свободной границы и распределение потенциала на ней для момента времени t0 + 2 и т.д.

Описанный алгоритм представляет собой явную схему Эйлера первого порядка точности, которая, как известно из теории числен ных методов, является условно устойчивой и требует наличия не которых дополнительных условий, связывающих шаг по времени с шагами по пространству (аналог условия Куранта в газовой дина мике [29]).

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

1) любая частица жидкости за временной шаг не может пере меститься на расстояние больше заданного;

2) узлы любого элемента не могут изменять ориентацию отно сительно друг друга (исключается самопересечение границы области).

Более подробно алгоритмы движения и выбора шага по време ни изложены в работах [14, 43].

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

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

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

Для использования в приложениях, составляющих основу па кета прикладных программ “AKORD”, был разработан оригиналь ный набор компонентов, позволяющих легко оперировать графиче скими объектами (рис. 7). Основными элементами этих компонен тов являются:

1. Графическое окно, которое выполняет следующие функции:

• преобразование логических (экранных) координат в физи ческие и обратно;

• произвольное масштабирование изображения и увеличение выделенных областей окна (функция “зуммирования”);

• автоматическая прорисовка осей координат с широкими функциями настройки;

• реализация быстрой прорисовки графики через страницы видеопамяти.

Глава 3. Математические методы и алгоритмы 2. Набор базовых графических объектов, на основе которых может быть реализован графический объект практически любой сложности. Сюда входят:

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

• вершина – объект, определяющий опорную точку для зада ния графического примитива. Количество опорных точек для задания объекта может быть жестко определено или варьироваться в зависимости от типа графического прими тива. Например, примитивы типа Line и Circle имеют две вершины, а примитив типа Arc - три, у примитива типа PolyLine количество вершин заранее не определено;

• граничная точка – объект, предназначенный для дискрети зации графического примитива по заданному закону. По умолчанию граничные точки совпадают с вершинами;

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

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

Глава 3. Математические методы и алгоритмы Рис. 7. Структура базовых графических объектов Базовые пространственные графические объекты Способы описания и работы с базовыми объектами Для однозначного задания пространственного объекта введем некоторые понятия. Пространственный объект задается, в первую очередь, последовательностью координат своих вершин. Кроме того, к последовательности вершин добавляются последовательно сти ребер пространственного объекта и его элементов.

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

Ребро представляется парой точек из исходной последователь ности точек, причем пары (m,l) и (l,m) (l и m - вершины из последо вательности вершин) считаются одним и тем же ребром. При опи Глава 3. Математические методы и алгоритмы сании элементов их ребра, так же как и вершины, задаются поряд ковыми номерами в последовательности номеров ребер, при этом каждое ребро упоминается не менее двух раз. В элементе ребра перечисляются в произвольном порядке.

Таким образом, пространственный объект определяется после довательностью номеров его вершин и ребер. Все дальнейшие преобразования выполняются над этими последовательностями.

Преобразования в пространстве Наилучшее восприятие формы пространственного объекта дос тигается при частом изменении его изображения в пространстве.

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

Для осуществления таких преобразований удобно пользоваться однородными координатами [70]. Точки в трехмерном пространст ве (x,y,z) представляются четырехмерным вектором (x,y,z,1) или (X,Y,Z,H).

Представление n-мерного вектора (n+1)-мерным назы вается однородным координатным воспроизведением [70]. При од нородном координатном воспроизведении конечные результаты в n-мерном пространстве получаются с помощью обратного преобра зования. То есть, трехмерный вектор (x,y,z) представляется четы рехмерным (hx, hy, hz, h). Разделив компоненты приведенного век тора на h, получим x=hx/h, y=hy/h, z=hz/h. Однородное координат ное воспроизведение позволяет осуществлять полный набор преоб разований в пространстве и, кроме того, дает возможность пред ставлять бесконечно удаленные точки векторами, состоящими из конечных компонент, например (1, 1, 1, 0). Необходимо отметить, что однородное координатное воспроизведение точек в трехмерном пространстве неединственно. Например, однородные координаты (12, 16, 8, 4), (6, 8, 4, 2) и (3, 4, 2, 1) представляют одну и туже точ ку (3, 4, 2).

Преобразование из однородных координат описывается соот ношением (X, Y, Z, H)=(x, y, z, 1)T, где T-матрица преобразования, которая может быть представлена в виде четырех отдельных частей Глава 3. Математические методы и алгоритмы 3 3 3 1 3 1 1.

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

a 0 0 b.

0 0 c 0 0 0 Помимо этого, диагональные элементы верхней левой подмат рицы 3х3 общей матрицы преобразования 4х4 позволяют осущест влять зеркальное отображение. Для отображения без изменения масштаба необходимо, чтобы определитель преобразования был равен -1. Так, например, при отображении относительно коорди натной плоскости XY изменяется только знак координаты z. Сле довательно, матрица преобразования для отображения относитель но плоскости XY имеет вид:

1 0.

0 0 0 00 Аналогично осуществляется отображение относительно других координатных плоскостей.

Недиагональные элементы подматрицы осуществляют сдвиг в трех измерениях, т.е.

Глава 3. Математические методы и алгоритмы 1 b c d 1 f z 1) = ( x + dy + gz bx + y + hz cx + fy + z 1) (x y g h 0 0 0.

Общая матрица вращения имеет вид n1 + (1 n1 )cos 2 n1n2 (1 cos ) + n3 sin n1n3 (1 cos ) n2 sin 2 n1n2 (1 cos ) n3 sin n2 + (1 n2 )cos n2 n3 (1 cos ) + n1 sin 2 n1n3 (1 cos ) + n2 sin n2 n3 (1 cos ) n1 sin n3 + (1 n3 )cos 0 0, где - угол вращения, а направление оси вращения представлено единичным вектором n = ( n1 n2 n3 ), (где n1, n2, n3 - косинусы уг лов, образованных осью вращения n с осями координат x, y и z, соответственно).

Матрицы вращения около осей координат получаются как ча стные случаи приведенной матрицы вращения около произвольной оси. Нижняя строка 1х3 общей матрицы преобразований произво дит перенос. В общем случае матрица переноса имеет вид 1 0 0 1.

0 0 l m n Правый столбец 3х1 осуществляет преобразования в перспек тиве. Последний скалярный элемент 1х1 выполняет общее измене ние масштаба. Полное преобразование, полученное путем воздей ствия на вектор матрицей положения 4х4 и нормализацией преоб разованного вектора, называется билинейным преобразованием [70] и обеспечивает выполнение комплекса операций сдвига, час Глава 3. Математические методы и алгоритмы тичного изменения масштаба, вращения, отображения, переноса, а также изменения масштаба изображения в целом.

Получение проекций Для получения проекции каждой точки P( x, y, z ), принадлежа щей визуализируемому объекту, предстоит вычислить точки изо бражения P( x, y, z ) на экране. Для этого необходимо преобразо вать координаты точки P из так называемых мировых координат ( x, y, z ) в экранные координаты ( x, y) ее центральной проекции.

Для получения экранных координат необходимо осуществить ви довое преобразование [2], при котором точка P остается на своем месте, но система мировых координат переходит в систему видо вых координат.

P P Рис. 8. Система видовых координат Для выполнения видовых преобразований должны быть заданы точки наблюдения и объекта (рис. 8). Выберем правую систему ко ординат. Кроме того, будем полагать, что экран находится между объектом и глазом E. Для точки P объекта прямая линия PE пере секает экран в точке P. Удобно, чтобы начало координат находи лось где-то вблизи центра объекта, так как объект наблюдается в направлении от Е к О.

Точку наблюдения удобно задать в сферических координатах E (,, ), то есть ее мировые координаты могут быть вычислены по формулам:

xe = sin cos, ye = sin sin, ze = cos.

Вектор ОЕ определяет направление наблюдения. Из точки на блюдения Е можно видеть точки объекта только внутри некоторого Глава 3. Математические методы и алгоритмы конуса, ось которого совпадает с линией ОЕ, а вершина - с точкой Е. При направлении взгляда из О в Е положительная полуось xe направляется вправо, а положительная полуось ye - вверх. Такое направление осей позволяет определить экранные координаты в тех же направлениях. Направление оси ze выбирается так, чтобы значение координат увеличивалось по мере удаления от точки на блюдения. При таком определении осей, которое представляется вполне удобным и логичным, система видовых координат будет левой. Тем не менее, такое несоответствие мировой и видовой сис тем координат не представляет никаких затруднений.

Видовое преобразование может быть записано в форме:

( xe ze 1) = ( xW 1) V, ye yW zW где V – матрица видового преобразования.

Видовое преобразование включает в себя следующую последо вательность действий.

1. Перенос системы координат, при котором точка Е становится но вым началом координат.

2. Поворот системы координат около оси Z на угол 2 в от рицательном направлении так, что ось Y совпадет по направле нию с горизонтальной составляющей отрезка ОЕ, а ось X будет располагаться перпендикулярно отрезку ОЕ.

3. Другой поворот около оси X на угол в положительном на правлении для того, чтобы ось z совпала с отрезком ОЕ.

4. Необходимо изменить ориентацию оси X, выполнив отображе ние относительно плоскости XZ.

Матрица V представляет собой произведение матриц перечис ленных линейных преобразований. После некоторого упрощения она будет иметь вид:

sin cos cos sin cos cos cos sin sin sin V =.

sin cos 0 0 0 0 Глава 3. Математические методы и алгоритмы После чего можно использовать видовые координаты xe и ye, отбросив координату ze. Ее вычисляют в том случае, когда необхо димо получить перспективные изображения и при последующей реализации алгоритма удаления невидимых линий. В этих случаях необходимо иметь полную картину видового преобразования.

Удаление невидимых линий и поверхностей Для стирания невидимых линий был выбран алгоритм Робертса [70], который представляет собой первое известное решение задачи об удалении невидимых линий. Алгоритм является векторным.

Это математически элегантный метод, работающий в объектном пространстве, то есть, используя нашу терминологию, в простран стве видовых координат. Алгоритм работает с выпуклыми телами.

Если тело невыпуклое, то оно перед применением этого алгоритма должно быть разбито на совокупность выпуклых тел.

Разрезание невыпуклых тел Ниже приводится алгоритм разбиения невыпуклых тел на вы пуклые части.

Процедура разрезания для каждой полигональной грани тела такова:

1. Перенести тело так, чтобы одна из вершин выбранной грани совпала с началом координат.

2. Повернуть тело вокруг начала координат так, чтобы одно из инцидентных ему ребер совпало с одной из осей координат, например, с осью x.

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

4. Проверить знаки координаты, которая перпендикулярна вы бранной грани (т.е. координаты z ), для всех вершин тела, не лежащих на этой грани.

5. Если все знаки совпадают или равны нулю, то тело является выпуклым относительно этой грани. В противном случае оно невыпукло;

разрезать тело плоскостью, несущей выбранную грань.

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

Глава 3. Математические методы и алгоритмы Метод Робертса для одного выпуклого тела Уравнение произвольной плоскости в трехмерном пространст ве имеет вид ax + by + cz + d = 0.

В матричной форме:

a b ( x y z 1) = 0, c d или ( x y z 1)( P )T = 0, где ( P ) = ( a b c d ) представляет собой плоскость. Поэтому любое выпуклое твердое тело можно выразить через матрицу, со стоящую из коэффициентов уравнений плоскостей, т.е.

a1 a2 … an b b... b (V ) = 1 2 n, c1 c2... cn d1 d 2... d n где каждый столбец содержит коэффициенты одной плоскости.

Напомним, что любая точка пространства представима в одно родных координатах вектором ( S ) = ( x y z 1). Более того, если точка ( S ) лежит на плоскости, то ( S )( P ) = 0. Если же ( S ) не ле T жит на плоскости, то знак этого скалярного произведения показы вает, по какую сторону расположена точка. В алгоритме Робертса предполагается, что точки, лежащие внутри тела дают положи тельное скалярное произведение.

Хотя уравнение плоскости содержит четыре неизвестных ко эффициента, его можно нормировать так, чтобы d = 1. Следова тельно, трех точек достаточно для определения этих коэффициен Глава 3. Математические методы и алгоритмы Подстановка координат этих точек, ( x1, y1, z1 ), тов.

( x2, y2, z2 ), ( x3, y3, z3 ) в нормированное уравнение дает:

z1 a x1 y x z2 b = 1.

y x z3 c y 3 Или в матричной форме:

( X )( C ) = ( D ). (3.27) Из решения этого уравнения определяются значения коэффициен тов уравнения плоскости:

( C ) = ( X )1 ( D ).



Pages:   || 2 | 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.