Механика сплошных сред (стр. 5 )

7.10.2. Краткая характеристика численных методов

7.10.2.1. Метод конечных разностей

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

Этот численный подход более всего развит в данное время и широко используется для решения как линейных, так и нелинейных уравнений гиперболического, эллиптического и параболического типа. Область интегрирования здесь разбивается на счетные ячейки с помощью некоторой, как правило, прямоугольной фиксированной сетки. Производные функции по всем направлениям заменяются конечными разностями с помощью тех или иных соотношений (приемы построения разностных уравнений весьма разнообразны). Причем, чаще всего используются так называемые неявные схемы. Тогда на каждом шаге приходится решать систему линейных алгебраических уравнений, содержащих иногда несколько тысяч неизвестных. Литература по этому направлению очень обширна (см., например [1-3]). Много внимания при этом уделяется исследованию свойств разностных уравнений (точность аппроксимации, условия устойчивости, диссипативные эффекты схем и т. п.).

7.10.2.2. Метод интегральных соотношений

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

7.10.2.3. Метод характеристик

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

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

7.10.2.4. Метод частиц в ячейках

Указанное направление численного моделирования сочетает в себе в определенных чертах преимущества лагранжева и эйлерова подходов. Область решения здесь разбивается неподвижной (эйлеровой) сеткой; однако сплошная среда трактуется дискретной моделью – рассматривается совокупность частиц фиксированной массы (лагранжева сетка частиц), которые и перемещаются через эйлерову сетку ячеек. Частицы служат для определения параметров самой жидкости (массы, энергии, скорости), в то время как эйлерова сетка используется для определения параметров поля течения (давления, плотности, температуры).

Метод частиц в ячейках позволяет исследовать сложные явления в динамике многокомпонентных сред, взаимодействия разрывов, поскольку частицы хорошо «следят» за свободными поверхностями и линиями раздела сред. Однако дискретный метод частиц обладает и рядом недостатков. Главный из них, лежащий в самой природе метода, состоит в том, что из-за дискретного представления сплошной среды (конечное число частиц в ячейке) методу присуща вычислительная неустойчивость (флуктуации). Затруднительно также получение информации для сильно разреженных областей, откуда практически уходят все частицы, и т. п.

7.10.2.5. Метод конечных элементов

В этом методе исходные уравнения и динамические краевые условия удовлетворяются только в некотором осредненном смысле для выбранного типичного конечного объема («элемента») среды. При этом аппроксимация различных полей проводится на конечном элементе локально и независимо от его положения в общей модели. Основная сфера приложения указанного подхода – это механика твердого деформированного тела. На основе данного метода построен известный программный комплекс ANSYS.

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

С другой стороны, методики этого типа имеют свою область применения и свои характерные трудности. По способу представления приближенного решения (которое обычно непрерывно или непрерывно с рядом производных) такие подходы, прежде всего, приспособлены для нахождения решения задач эллиптического и параболического типов. При решении гиперболических задач методы конечных элементов нельзя считать достаточно эффективными. Основная причина заключается в том, что здесь полностью отсутствует использование такого фундаментального свойства гиперболических задач, как конечность области влияния. Это приводит к неестественному «завязыванию» всех узлов расчетной области, следствием чего являются неоправданно высокие (для задач гиперболического типа) требования к объему используемой памяти ЭВМ.

7.10.2.6. Метод дискретных вихрей

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

7.10.2.7. Статистические методы

Быстрое развитие вычислительной техники стимулировало разработку численных методов статистического моделирования (методы Монте-Карло) широкого класса задач механики жидкости, физики, биологии, химии. Этот класс задача условно можно разделить на два вида:

1. Задачи со стохастической природой. Для данных задач метод Моне-Карло используется для прямого моделирования естественной вероятностной модели. При этом точная динамика заменяется стохастичным многомерным процессом;

2. Детерминированные задачи. Указанные задачи описываются вполне определенными уравнениями. Здесь искусственно строится вероятностный процесс, который численно моделируется методом Монте-Карло на ЭВМ, что позволяет получить формальное решение в виде статистических оценок. При этом необходимо показать адекватность построенного вероятностного процесса рассматриваемому кинетическому уравнению.

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

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

7.10.3. Основы численных методов

7.10.3.1. Задача интерполирования

При вычислениях оперируют с сеточными функциями, т. е. функциями, заданными на дискретной совокупности точек – узлов сетки [1-3]. Если нужно знать значения при , не совпадающих с узлами, то поступают следующим образом. Строят некоторую достаточно простую функцию , которая совпадет с в узлах . В промежуточных значениях функция приближенно представляет функцию . Эту функцию называют интерполирующей, а задачу ее отыскания – интерполированием.

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

7.10.3.2. Интерполяционный многочлен Лагранжа

Рассмотрим простейший и самый распространенный случай, когда - многочлен. Требуется построить многочлен степени , который в заданных точках (узлах интерполирования) принимает заданные значения . Легко доказать, что такой многочлен только один. При помощи простых выкладок можно построить интерполяционный многочлен вида [5]

, (7.10.1)

где

. (7.10.2.)

Интерполяционный полином (7.10.1) называют интерполяционным полиномом Лагранжа.

7.10.3.3. Погрешность интерполирования

Для получения информации о погрешности интерполирования необходимо оценить разность , считая функцию достаточно гладкой. Можно получить [5]:

(7.10.3)

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

Формула (7.10.3) позволяет оценить погрешность полиномиальной интерполяции, если известна оценка для производной . Полагая

, (7.10.4)

имеем

. (7.10.5)

7.10.4. Вычисление интегралов

7.10.4.1. Квадратурные формулы Ньютона-Котеса

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

, (7.10.6)

где оценивается по формуле (7.10.3). Тогда имеем

, (7.10.7)

где - погрешность интегрирования. Из (7.10.5) следует, что

. (7.10.8)

Подставляя (7.10.1) в (7.10.7), получим

, (7.10.9)

где

(7.10.10)

- коэффициенты квадратурной формулы. Формулы (7.10.9) называют формулами Ньютона-Котеса. Рассмотрим частные случаи.

7.10.4.2. Формула трапеций

Пусть , тогда , , Коэффициенты легко вычислить: Таким образом

. (7.10.11)

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

. (7.10.12)

7.10.4.3. Формула Симпсона

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

. (7.10.13)

Формулу (7.10.13) называют формулой Симпсона. Разделив отрезок на частей и применив к интервалам формулу (7.10.13), получаем

. (7.10.14)

7.10.5. Численное дифференцирование

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

, (7.10.15)

. (7.10.16)

Рассмотрим случай при равноотстоящих узлах , . Легко показать, что в этом случае выражения для вычисления производной будут выглядеть так:

, (7.10.17)

(7.10.18)

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

Формулы для производных при имеют вид:

,

(7.10.19)

Порядок аппроксимации этих формул равен двум.

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

7.11. Применение метода потоков в механике сплошных идеальных сред

7.11.1. Общие замечания

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

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

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

Не случайно, видимо, при больших значениях Re решение полных уравнений Навье – Стокса с сохранением влияния членов молекулярной вязкости весьма затруднительно. Основная трудность, возникающая при их применении, состоит в достаточно точном «разрешении» структуры потока при не слишком малых размерах расчетной сетки (шаг сетки должен быть таким, чтобы погрешность аппроксимации конвективных членов была бы много меньше разностных представлений вязкостных членов). Эта трудность может быть частично преодолена применением сгущающихся в нужных местах сеток и схем повышенной точности.

При расчете таких моделей на реальных («грубых») сетках формальное решение может быть получено и для больших значений чисел Рейнольдса. Однако, такое решение может, вообще говоря, не соответствовать уравнениям Навье – Стокса, так как молекулярные эффекты здесь могут «забиваться» схемной (эффективной) вязкостью, обеспечивающей вычислительную устойчивость решения в целом. Данный подход (с приближенным механизмом диссипации энергии) можно использовать лишь для задач, где влияние вязкости незначительно и течение автомодельно по Re. Таким образом, представляется важным, чтобы алгоритм расчета вязких течений позволял осуществить предельный переход к моделям идеального газа, когда кинематическая вязкость .

Далее описывается один из численных подходов к изучению свойств течений сжимаемого газа – метод потоков. Отличительная черта указанного подхода заключается в численном решении с помощью консервативных разностных схем общих нестационарных уравнений, записанных в виде законов сохранения в интегральной форме для идеального и вязкого теплопроводного сжимаемого газа. Такой подход позволяет, по существу, избежать трудностей, связанных с аппроксимацией старших производных уравнений Навье – Стокса. Разделение вектора плотности потока на конвективную и вязкую составляющие дает возможность использовать данную методику для расчета движений, как вязкой, так и идеальной среды.

7.11.2. Описание метода потоков

В данном разделе будут изложены общие принципы построения конечно-разностных схем метода потоков, разработанного и , [6].

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

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

Зная значения и Е можно вычислить средние для данной ячейки величины плотностей распределения перечисленных количеств:

, , . (7.11.1)

где - массовая плотность газа; - i – компонента плотности распределения импульса; - плотность распределения полной энергии.

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

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

, . (7.11.2)

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

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

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

, (7.11.3)

где - вектор плотности потока величины , - единичный вектор внешней к нормали. Векторы плотностей потоков определяются на поверхности переменными газодинамического поля и их производными и плотностями распределения. Уравнения (7.11.3) справедливы для произвольного объема, и естественно требовать их выполнения для минимального объема-ячейки разностной сетки.

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

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

При численном решении интегралы в (7.11.3) вычисляются по квадратурной формуле

, (7.11.4)

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

Уравнения (7.11.4), выписанные для каждой аддитивной характеристики среды во всех элементарных объемах (ячейках), составляют систему нелинейных алгебраических уравнений относительно переменных газодинамического поля в одной точке элементарного объема (ячейки) . Система уравнений (7.11.4) вместе с соотношениями (7.11.1), (7.11.2), дополненная уравнением состояния и граничными условиями, является замкнутой.

Остановимся на отличительных чертах рассматриваемой численной модели. Все переменные в (7.11.3) естественно разделяются на две группы: первую – переменные газодинамического поля (составляющие вектора скорости , давление , температура ) и вторую – переменные плотности распределения количеств - потоки, для которых сформулированы законы сохранения (7.11.3).

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

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

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

Переменные же газодинамического поля и их производные (в тензоре вязких напряжений и законе теплопроводности) определяются по симметричным формулам. На аппроксимацию переменных газодинамического поля накладывают ограничения и другие законы механики.