Решение уравнения переноса методом конечных элементов

Численное решение методом конечных элементов задач диффузионного и конвективного переноса в сильно гетерогенных пористых средах Текст научной статьи по специальности « Математика»

Аннотация научной статьи по математике, автор научной работы — Васильева Мария Васильевна, Васильев Василий Иванович, Тимофеева Татьяна Семеновна

В работе рассматривается конечно-элементная аппроксимация уравнения конвективного и диффузионного переноса. Рассмотрены различные методики стабилизации конечно-элементной аппроксимации по пространственным переменным, такие как противопотоковая аппроксимация конвективного слагаемого посредством введения искусственной диффузии (AD, artificial diffusion), метод SUPG (streamline upwind Petrov-Galerkin), которые используются для стабилизации классического метода Галеркина . В работе рассмотрена также аппроксимация уравнения переноса с использованием разрывного метода Галеркина , позволяющая строить противопотоковые схемы аппроксимации. Представлены результаты численного сравнения рассматриваемых схем на примере задачи конвективного и диффузионного переноса в пористых средах. Рассмотрены задачи фильтрации с сильно гетерогенными коэффициентами, которые ведут к большим перепадам давления (большим скоростям фильтрации , градиентам давления).

Похожие темы научных работ по математике , автор научной работы — Васильева Мария Васильевна, Васильев Василий Иванович, Тимофеева Татьяна Семеновна

The finite element approximation of the convective and diffusive transport equation has been considered. Different methods for stabilization of the finite element approximation have been discussed: upwind approximation of the convective term using artificial diffusion (AD) and streamline upwind Petrov-Galerkin (SUPG) method, both used for stabilization of the classic Galerkin method . Another approach to approximation of the transport equation related to the discontinuous Galerkin method (DG) has been investigated. This method also allows to approximate the convective term using upwind schemes. The results of the numerical comparison of the considered schemes for the convective and diffusive transport problems in a porous media have been presented. The flow and transport in a highly contrast heterogeneous porous media that lead to the significant pressure gradients and, therefore, high velocities have been considered as test problems.

Текст научной работы на тему «Численное решение методом конечных элементов задач диффузионного и конвективного переноса в сильно гетерогенных пористых средах»

2016, Т. 158, кн. 2 С. 243-261

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

ISSN 1815-6088 (Print) ISSN 2500-2198 (Online)

ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ ЗАДАЧ ДИФФУЗИОННОГО И КОНВЕКТИВНОГО ПЕРЕНОСА В СИЛЬНО ГЕТЕРОГЕННЫХ ПОРИСТЫХ СРЕДАХ

М.В. Васильева1,2, В.И. Васильев1, Т.С. Тимофеева1

1 Северо-Восточный федеральный университет имени М.К. Аммосова, г. Якутск, 677000, Россия 2Институт вычислительной математики и математической геофизики СО РАН,

г. Новосибирск, 630090, Россия Аннотация

В работе рассматривается конечно-элементная аппроксимация уравнения конвективного и диффузионного переноса. Рассмотрены различные методики стабилизации конечно-элементной аппроксимации по пространственным переменным, такие как проти-вопотоковая аппроксимация конвективного слагаемого посредством введения искусственной диффузии (AD, artificial diffusion), метод SUPG (streamline upwind Petrov-Galerkin), которые используются для стабилизации классического метода Галеркина. В работе рассмотрена также аппроксимация уравнения переноса с использованием разрывного метода Галеркина, позволяющая строить противопотоковые схемы аппроксимации. Представлены результаты численного сравнения рассматриваемых схем на примере задачи конвективного и диффузионного переноса в пористых средах. Рассмотрены задачи фильтрации с сильно гетерогенными коэффициентами, которые ведут к большим перепадам давления (большим скоростям фильтрации, градиентам давления).

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

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

Изучаемые процессы описываются системой уравнений для скорости и давления, а также уравнением конвективного и диффузионного переноса для концентрации примеси [3, 6, 8, 9]. Рассмотрим вычислительные алгоритмы решения поставленной задачи, основанные на методе конечных элементов [10]. Для аппроксимации системы уравнений относительно скорости и давления используется смешанный метод конечных элементов [11—13].

Основные сложности при численном решении связаны с уравнением переноса [14, 15]. Классический метод Галеркина, когда тестовые и триальные базисные функции совпадают, в случае задач с доминирующим конвективным слагаемым приводит к численной неустойчивости — возникновению осцилляций в решении задачи, поэтому для таких задач обычно применяют специальные методы стабилизации [16, 17]. Среди классических методов стабилизации можно выделить такие методы, как противопотоковая аппроксимация конвективного слагаемого посредством введения искусственной диффузии (AD, artificial diffusion), метод SUPG (streamline upwind Petrov-Galerkin). Такие методики стабилизации подавляют нефизичные осцилляции в приближенном решении задачи. Основным недостатком перечисленных методов является появление дополнительных стабилизирующих членов в уравнении, которые приводят к усложнению структуры матрицы, кроме того, метод чувствителен к выбору параметров стабилизации, неправильный выбор которых приводит либо к потери устойчивости, либо к ухудшению точности [18]. Другим известным методом для аппроксимации уравнений конвективного и диффузионного переноса является разрывный метод Галеркина, который позволяет строить противопотоковые схемы аппроксимации для конвективного слагаемого уравнения 19.

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

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

Далее, в третьем разделе рассмотрен разрывный метод Галеркина для аппроксимации уравнения переноса, позволяющий строить противопотоковые аппроксимации конвективного слагаемого в уравнении переноса.

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

1. Постановка задачи и конечно-элементная аппроксимация

Уравнение переноса конвективного и диффузионного переноса имеет вид дс

——+ u • grad с — div(a grad с) = f, x G Q, T > 0, (1)

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

Уравнение (1) дополняется начальным условием

c(x, 0) = с0, x G Q, (2)

и граничными условиями

с = ci, x G Г1, с = С2, x G Г2,

— а— =0, x G OQ/Г /Г2, On

где n — вектор внешней нормали.

Для описания течения в пористой среде будем использовать следующую систему уравнений:

u = —k gradp, x G Q

Система уравнений (4) дополняется соответствующими граничными условиями:

p = pi, x G Г1, p = P2, x G Г2,

u = 0, x G 0Q/r1 /Г2,

Далее будем предполагать, что Q — прямоугольная область, а Г1 и Г2 — левая и правая границы соответственно.

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

Проведем дискретизацию линейного уравнения конвективного и диффузионного переноса (1) с помощью классического метода Галеркина, в котором тестовые и триальные функции совпадают. Пусть с и z G Q(Q), вариационную постановку запишем следующим образом:

-zdx + (aVen+1, Vz) dx + u •Ven+1zdx = 0 Vz G Q, (6)

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

Отметим, что для уравнения (6) сеточное число Пекле определяется следующим образом: Peh = (umaxh)/(2amax), где amax — верхняя граница коэффициента диффузии, umax — максимальная скорость фильтрации и h — характерный шаг сетки. При малых числах Пекле Peh 1 представленная аппроксимация может приводит к появлению осцилляций в решении, и поэтому необходимо использовать специальные схемы аппроксимации со стабилизацией.

Для дискретизации системы уравнений (4) применим смешанный метод конечных элементов и получим следующую вариационную формулировку задачи: найти (и,р) € V(О) х Ь(О) такие, что

где (V, д) € V (О) х Ь(О), с! = 2, 3 .В качестве базисных функций для поля скоростей будем использовать элементы Равьяра-Тома (ИТ, Ил^ай-ТЬотаз) [12] и кусочно-постоянные функции для давления, удовлетворяющие ЬББ (ш^ир) условию [11].

2. Аппроксимация уравнения переноса со стабилизацией

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

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

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

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

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

ь Уз — Уз-1 — а Уз+1 — 2Уз + Уз-1 = ь Уз+1 — Уз-1 — ( + Уз+1 — 2уз + Уз-1

h h2 2h V 2 J h2

где b и a — коэффициенты конвекции и диффузии, а h — шаг сетки. Таким образом, аналогичная аппроксимация может быть получена при использовании центрально-разностной аппроксимации уравнения bux — (a + bh/2) uxx = 0.

В методе конечных элементов аналогичная технология может быть также использована для достижения противопотоковых аппроксимаций [16, 17]. Тогда вместо уравнения (1) будем рассматривать модифицированное уравнение c дополнительным слагаемым (искусственная диффузия, AD)

— + u ■ Ус — div (âgrad с) = 0, (8)

где â = (a + \u\h/2) есть модифицированный коэффициент диффузии с дополнительным слагаемым (искусственная диффузия), который пропорционален шагу по сетке и максимальной скорости фильтрации, \u\ = \J(u, u).

Таким образом, аппроксимацию уравнения (8) с использованием метода конечных элементов запишем в виде: найти cn+1 G Q такую, что

-zdx + (â Vcn+1, Vz) dx + u ■Vcn+1zdx = 0 V z G Q, (9)

Отметим, что подобная аппроксимация имеет первый порядок в уравнении переноса, вместо искомого уравнения (1) мы будем решать модифицированное уравнение (8).

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

В рассмотренном ранее методе добавляется искусственная диффузия, равная aAD = \u\h/2 (изотропный коэффициент диффузии). Пусть теперь будем добавлять искусственную диффузию только в направлении течения e = u/\u\, то есть ‘a = aAD ee = a ad (uu)/\u\2 (анизотропный коэффициент диффузии). Тогда дополнительное слагаемое можно записать следующим образом:

I (aVcn+1, Vz) dx = j(aAD U Vcn+1, Vz) dx = nn

( (uVcn+1,0AD uVAdx = ¡u ■Vcn+lhvl ■ Vz dx. (10) J V \u\2 J J 2\u\ V ;

Таким образом, аппроксимацию уравнения переноса запишем в виде

-zdx + (аЧсп+1,Чг) dx+

+ f u ■Vcn+1hu ■Vzdx + f u ■Vcn+1zdx = 0, (11)

-zdx + (aVcn+1, Vz) dx + u ■ Vcn+1Z dx = 0 (12)

h , где z = (z + ——гu ■ Vz).

Вариационная постановка в виде (12) называется методом SU (streamline upwind) [16, 17]. Основным минусом такой постановки является то, что уравнение уже не соответствует исходному дифференциальному уравнению, поэтому вместо уравнения (12) в методе SUPG используется следующая вариационная постановка: найти cn+1 е Q

/ cn+1 — cn Г Г _ -Zdx + (aVcn+1, Vz) dx + u ■Vcn+1Zdx = 0 V Z е Q, (13)

Таким образом, поскольку триальное пространство Q и тестовое пространство Q не совпадают, этот метод является методом Петрова -Галеркина (PG, Petrov-Galerkin). Полученная аппроксимация уже не может быть интерпретирована как схема с искусственной диффузией, поскольку само уравнение остается прежним, а меняется только тестовая функция z .

3. Аппроксимация уравнения переноса с использованием разрывного метода Галеркина

Разрывный метод Галеркина был впервые предложен в [24] для решения гиперболического уравнения, с тех пор метод активно развивается. В [22, 23] метод LDG (local discontinuous Galerkin) рассматривается для решения задач конвекции-диффузии. Метод в различных его вариациях активно разрабатывается и для решения эллиптических уравнений [19, 21, 25], для задач теории упругости [19, 26], а также для задач гидродинамики 28. Метод IPDG (interior penalty discontinious Galerkin) рассмотрен в работах [19, 20, 31], в том числе и для решения задач многофазной фильтрации. Отметим что метод IPDG, который мы и будем использовать в настоящей работе, имеет несколько вариаций: NIPG (nonsymmetric IPDG), SIPDG (symmetric IPDG), IIPDG (incomplete IPDG) [21, 32, 33].

В статье [21] рассмотрено обобщение разрывного метода Галеркина, где предложены различные варианты метода для эллиптического оператора. Разрывный метод Галеркина также применяется для построения аппроксимаций по времени и рассматривается как некоторый общий случай, объединяющий такие классические аппроксимации, как явная, неявная схемы (backward and forward Euler) для аппроксимации по времени [34], методы высокого порядка Рунге-Кутта [23]. Метод интересен для рассмотрения и является более общей методологией аппроксимации не только по пространству, но и по времени. В настоящее время метод очень активно развивается, и на его основе создается множество других современных, интересных методов, такие как DPG, HDG и др. 37. В [39, 40] разработана абстрактная теория схем разрывного метода Галеркина в смешанной формулировке для квазилинейных эллиптических уравнений второго порядка.

При сравнении со стандартным методом Галеркина (CG, continious Galerkin) следует отметить, что разрывный метод Галеркина является более современным и в настоящее время активно развивается и применяется для решения задач различного типа. А классический метод Галеркина имеет более чем 60-летнюю историю, и о нем написано большое количество работ и книг, что, конечно, является преимуществом, но все же представляет научный интерес использование и рассмотрение современных методов, активно разрабатываемых в научной среде, которые основаны

на разрывном методе Галеркина (DG). Одним из достоинств метода DG является локальная консервативность аппроксимации, в то время как метод CG обладает только глобальной консервативностью. Свойство локальной консервативности для задач переноса является критическим, например, для решения задач теории фильтрации. Среди недостатков метода DG следует отметить большую размерность задачи — количество неизвестных получаемой системы уравнений. А именно, в методе DG при аппроксимации с использованием линейных базисных функций (полиномов первой степени) размерность задачи (количество степеней свободы, DOF) будет равняться произведению количества элементов на количество узлов в элементе, а в методе CG размерность задачи равна количеству узлов расчетной сетки, что существенно меньше. Например, для структурированной треугольной сетки 5 на 5, то есть содержащей 50 треугольных элементов и 36 узлов, размерность задачи при аппроксимации с использованием метода CG будет N = 36, а для метода DG (IPDG) — N = 50 ■ 3 = 150. И, следовательно, время счета может быть существенно больше.

Рассмотрим аппроксимацию уравнения (1) с использованием разрывного метода Галеркина, в частности SIPDG (symmetric interior penalty discontinious Galerkin). Разрывный метод Галеркина позволяет также строить противопотоковые аппроксимации конвективного слагаемого. Пусть Th — триангуляция области О и Г^ -множество всех внутренних граней. Пусть e есть грань между двумя элементами сетки K1 и K2 (треугольниками или тетраэдрами), тогда определим функцию скачка и среднего для функции u следующим образом:

М = uKl + , [u] = u|K ■ п\кг — u\k2 ■ п\к2. (14)

Вариационная постановка для задачи (1) с использованием разрывного метода Галеркина IPDG записывается следующим образом: найти c е QDG такую, что

J2 —zdx + aDG(c, z) = Idg(v) Vz е QDG, (15)

aDG(u, v) = / aVcVzdx — ^^ [z] ds — ^^ [c] ds+

кетh K ee_rh e ee_rh <

+ Y / a[c][z] ds + / u ■ necup[z] ds + / u ■ neczds, (16) eerh < eeThJe eer2 <

Idg(v) = / fzdx + u ■ nec1 z ds (17)

и сир = сК1 при и ■ пе > 0 и сир = сК2 при и ■ пе Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

aM.K. Ammosov North-Eastern Federal University, Yakutsk, 677000 Russia bInstitute of Computational Mathematics and Mathematical Geophysics, Siberian Branch, Russian Academy of Sciences, Novosibirsk, 630090 Russia E-mail: *vasilyevadotmdotv@gmail.com, **vasvasil@mail.ru, ***timofeevatc52@mail.ru

Received March 23, 2016 Abstract

The finite element approximation of the convective and diffusive transport equation has been considered. Different methods for stabilization of the finite element approximation have been discussed: upwind approximation of the convective term using artificial diffusion (AD) and streamline upwind Petrov-Galerkin (SUPG) method, both used for stabilization of the classic Galerkin method. Another approach to approximation of the transport equation related to the discontinuous Galerkin method (DG) has been investigated. This method also allows to approximate the convective term using upwind schemes. The results of the numerical comparison of the considered schemes for the convective and diffusive transport problems in a porous media have been presented. The flow and transport in a highly contrast heterogeneous porous media that lead to the significant pressure gradients and, therefore, high velocities have been considered as test problems.

Keywords: convection-diffusion equation, filtration, heterogeneous porous media, finite element method, numerical stabilization, classic Galerkin method, discontinuous Galerkin method

Acknowledgments. This study was supported in part by the Russian Science Foundation (project no. 15-11-10024) and the Russian Foundation for Basic Research (project no. 15-3120856).

Fig. 1. Heterogeneous permeability coefficient: on the left -variant 1, on the right — variant 2.

Fig. 2. Filtration velocities ux (from above) and uy (from below): on the left — variant 1, on the right — variant 2.

Fig. 3. Solution of the convective and diffusive transport problem at a = 0.001 (Task 1) for variant 1 (on the left) and variant 2 (on the right) using the classic Galerkin method (CG) with stabilization by artificial diffusion (AD), streamline upwind Petrov-Galerkin (SUPG) method, and discontinuous Galerkin method (DG) (from up to down) at the finite moment of time.

Fig. 4. Isolines for c = 0.9 of solution of the convective and diffusive transport problem (Problem 1) for variant 1 (on the left) and variant 2 (on the right) for different diffusive transport coefficients (from up to down) using various space approximation methods at the finite

moment of time. Classic Galerkin method (CG) — black line; with stabilization by artificial diffusion (AD) — red line; SUPG method — green line; discontinuous Galerkin method (DG) -blue line.

Fig. 5. Solution of the convective and diffusive transport problem (Problem 1 ) along the line y = 0.18 for variant 1 (from above) and variant 2 (from bottom) for different diffusive transport coefficients (from right to left) at the finite moment of time. Classic Galerkin method (CG) — black line; with stabilization by artificial diffusion (AD) — red line; SUPG method -green line; discontinuous Galerkin method (DG) — blue line.

Fig. 6. Solution of the convective transport problem (Problem 2) for variant 1 (on the left) and variant 2 (on the right) at different moments of time tm , m = 20, 60 and 100 (from up to down) using the SUPG method.

Fig. 7. Solution of the convective transport problem (Problem 2) at the finite moment of time for variant 1 (on the left) and variant 2 (on the right) using the classic Galerkin method (CG), as well as stabilization by artificial diffusion (AD), SUPG method, and discontinuous Galerkin method (DG) (from up to down).

1. Aziz K., Settari A. Petroleum Reservoir Simulation. Appl. Sci. Publ. Ltd., 1979. 497 p.

2. Bear J. Dynamics of Fluids in Porous Media. N. Y., Dover Publ., Inc., 1972. 764 p.

3. Chen Z., Huan G., Ma Y. Computational Methods for Multiphase Flows in Porous Media. SIAM, 2006. 561 p.

4. Vabishchevich P.N., Vasil’eva M. Iterative solution of the pressure problem for the multiphase filtration. Math. Model. Anal., 2012, vol. 17, no. 4, pp. 532-549. doi: 10.3846/13926292.2012.706655.

5. Christie M.A., Blunt M.J. Tenth SPE comparative solution project: a comparison of upscaling techniques. SPE Reservoir Eval. Eng., 2001, vol. 4, no. 4, pp. 308-317.

6. Chung E.T., Efendiev Y., Lee C.S. Mixed generalized multiscale finite element methods and applications. Multiscale Model. Simul., 2015, vol. 13, no. 1, pp.338-366. doi: 10.1137/140970574.

7. Chung E.T., Efendiev Y., Leung W.T., Ren J. Multiscale simulations for coupled flow and transport using the generalized multiscale finite element method, Computation, vol. 3, no. 4, pp.670-686. doi: 10.3390/computation3040670.

8. Aarnes J.E., Gimse T., Lie K.A. Geometric Modelling, Numerical Simulation, and Optimization. An introduction to the numerics of flow in porous media using Matlab. Berlin, Heidelberg, Springer, 2007, pp.265-306.

9. Jenny P., Lee S.H., Tchelepi H.A. Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media. J. Comput. Phys., vol. 217, no. 2, pp. 627-641. doi: 10.1016/j.jcp.2006.01.028.

10. Brenner S., Scott R. The Mathematical Theory of Finite Element Methods. Springer Science & Business Media, 2007. 400 p.

11. Brezzi F., Fortin M. Mixed and Hybrid Finite Element Methods. Berlin: Springer, 1991. 350 p.

12. Raviart P.A., Thomas J.M. Mathematical Aspects of Finite Element Methods. A mixed finite element method for 2-nd order elliptic problems. Berlin, Heidelberg,Springer, 1977, pp. 292-315.

13. Carstensen C. A posteriori error estimate for the mixed finite element method. Math. Comput., Am. Math. Soc., 1997, vol. 66, no. 218, pp. 465-476.

14. Vabishchevich P.N., Vasil’eva M. V. Explicit-implicit schemes for convection-diffusion-reaction problems. Numer. Anal. Appl., 2012, vol. 5, no. 4, pp. 297-306. doi: 10.1134/S1995423912040027.

15. Afanas’eva N.M., Vabishchevich P.N., Vasil’eva M.V. Unconditionally stable schemes for convection-diffusion problems. Russ. Math., 2013, vol. 57, no. 3, pp. 1-11. doi: 10.3103/S1066369X13030018.

16. Donea J., Huerta A. Finite Element Methods for Flow Problems. Chichester, John Wiley & Sons Ltd., 2003. 362 p.

17. Brooks A.N. A Petrov-Galerkin finite element formulation for convection dominated flows. Ph.D. Diss., California Institute of Technology, 1981. 126 p.

18. Ol’shanskii M.A. An analysis of the multigrid method for the convection-diffusion equations with the Drichlet boundary conditions. Comput. Math. Math. Phys., 2004, vol. 44, no. 8, pp. 13741403.

19. Riviere B. Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation. SIAM, 2008. 212 p.

20. Riviere B., Wheeler M.F. Discontinuous Galerkin methods for flow and transport problems in porous media. Commun. Numer. Methods Eng., 2002, vol. 18, no. 1. pp. 63-68.

21. Arnold D.N., Brezzi F., Cockburn B., and Marini L.D. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 2002, vol. 39, no. 5, pp. 1749-1779.

22. Cockburn B., Shu C.W. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal., 1998, vol. 35, no. 6, pp. 2440-2463.

23. Cockburn B., Shu C.W. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. J. Sci. Comput., 2001, vol. 16, no. 3. pp. 173-261.

24. Reed W.H., Hill T.R. Triangular Mesh Methods for the Neutron Transport Equation. Report LA-UR-73-479. Los Alamos, Los Alamos Sci. Lab., 1973. 23 p.

25. Brezzi F., Manzini G., Marini D., Pietra P., Russo A. Discontinuous Galerkin approximations for elliptic problems. Numer. Methods Partial Differ. Equations, 2000, vol. 16, no. 4, pp. 365-378.

26. Riviere B., Shaw S., Wheeler M.F., Whiteman J.R. Discontinuous Galerkin finite element methods for linear elasticity and quasistatic linear viscoelasticity. Numer. Math., 2003, vol. 95, no. 2, pp. 347-376.

27. Girault V., Riviere B., Wheeler M.F. A discontinuous Galerkin method with non-overlapping domain decomposition for the Stokes and Navier-Stokes problems. Math. Comput., 2005, vol. 74, no. 249, pp. 53-84.

28. Girault V., Riviere B. DG approximation of coupled Navier-Stokes and Darcy equations by Beaver-Joseph-Saffman interface condition. SIAM J. Numer. Anal., 2009, vol. 47, no. 3, pp. 2052-2089.

29. Cockburn B., Kanschat G., Schotzau D. The local discontinuous Galerkin method for the Oseen equations. Math. Comput., 2004, vol. 73, no. 246, pp. 569-593.

30. Cockburn B., Kanschat G., Schotzau D. A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations. J. Sci. Comput., 2007, vol. 31, nos. 1-2, pp. 61-73.

31. Sun S., Wheeler M.F. Discontinuous Galerkin methods for coupled flow and reactive transport problems. Appl. Numer. Math., 2005, vol. 52, no. 2, pp. 273-298.

32. Babuska I., Baumann C.E., Oden J.T. A discontinuous hp finite element method for diffusion problems: 1-D analysis. Comput. Math. Appl., 1999, vol. 37, no. 9, pp. 103-122.

33. Wheeler M.F. An elliptic collocation-finite element method with interior penalties. SIAM J. Numer. Anal., 1978, vol. 15, no. 1, pp. 152-161.

34. Gottlieb S., Wei G.W., Zhao S. A Unified Discontinuous Galerkin Framework for Time Integration. Preprint. 2010. 42 p. Available at: http://fodava.gatech.edu/files/reports/F0DAVA-10-34.pdf.

35. Kirby R.M., Sherwin S.J., Cockburn B. To CG or to HDG: a comparative study. J. Sci. Comput., 2012, vol. 51, no. 1, pp. 183-212. doi: 10.1007/s10915-011-9501-7.

36. Nguyen N.C., Peraire J., Cockburn B. An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations. J. Comput. Phys., 2009, vol. 228, no. 9, pp. 3232-3254. doi: 10.1016/j.jcp.2009.01.030.

37. Demkowicz L., Gopalakrishnan J. Analysis of the DPG method for the Poisson equation. SIAM J. Numer. Anal. 2011, vol. 49, no. 5, pp. 1788-1809. doi: 10.1137/100809799.

38. Demkowicz L., Demkowicz L., Gopalakrishnan J. A class of discontinuous Petrov-Galerkin methods. Part I: The transport equation. Comput. Methods Appl. Mech. Eng., 2010, vol. 199, no. 23, pp. 1558-1572.

39. Dautov R.Z., Fedotov E.M. Discontinuous mixed penalty-free Galerkin method for second-order quasilinear elliptic equations. Comput. Math. Math. Phys., 2013, vol. 53, no. 11, pp. 1614-1625.

40. Dautov R.Z., Fedotov E.M. Abstract theory of hybridizable discontinuous Galerkin methods for second-order quasilinear elliptic problems. Comput. Math. Math. Phys., 2014, vol. 54, no. 3, pp. 474-490.

41. Anders Logg, Kent-Andre Mardal, Garth N. Wells Automated Solution of Differential Equations by the Finite Element Method. The FEniCS Book. 2011. 723 p.

42. Software GMSH. Available at: http://geuz.org/gmsh/.

43. Software package PARAVIEW. Available at: http://www.paraview.org/.

Для цитирования: Васильева М.В., Васильев В.И., Тимофеева Т.С. Численное решение методом конечных элементов задач диффузионного и конвективного переноса \ в сильно гетерогенных пористых средах // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. — 2016. — Т. 158, кн. 2. — С. 243-261.

Лекция 10. Метод конечных элементов. Плоская задача.

10.1. Содержание метода.Метод конечных элементов (МКЭ) представляет

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

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

МКЭ предполагает иной подход. Рассматривается элемент конечных размеров (КЭ), за счет чего осуществляется переход от сплошной системы с бесконечным числом степеней свободы, к системе с конечным числом степеней свободы.

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

Поставив в соответствие каждому узловому усилию узловое перемещение, представим сплошной диск набором конечных элементов (КЭ), взаимодействующих между собой в конечном числе узловых точек.

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

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

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

Алгоритм решения задач по МКЭ содержит следующие этапы:

1. Дискретизация — разбиение заданной области на КЭ; нумерация узлов и КЭ.

2. Аппроксимация перемещений узлов в КЭ.

3. Построение матриц жесткости (МЖ) конечных элементов.

4. Построение глобальной матрицы жесткости общей системы уравнений.

5. Сведение нагрузок и воздействий, приложенных к КЭ, к узловым силам, учет условий закрепления.

6. Решение общей системы уравнений.

7. Определение напряжений и (при необходимости) деформаций в КЭ.

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

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

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

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

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

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

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

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

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

Ширина полосы B вычисляется по формуле

где R – максимальная по элементам величина наибольшей разности между номерами узлов в отдельном элементе; Q — число неизвестных (число степеней свободы в каждом узле). Минимизация величины В связана с минимизацией R, что может быть осуществлено последовательной нумерацией узлов при движении в направлении наименьшего размера тела против часовой стрелки. Правильная нумерация узлов экономит машинную память более чем на 60%, хотя не влияет на вычислительные аспекты задачи.

Построение градиентной матрицы.Для плоского КЭ «e» в форме треугольника с вершинами i, j, k с узловыми координатами соответственно ( , ), ( , ), ( , ), перемещения каждого узла в направлении осей x и y имеют две компоненты u и v.

= , (10.2)

а для трех узлов i, j, k шесть компонент образуют вектор перемещений внутри КЭ «e»:

= = (10.3)

Простейшими интерполяционными функциями являются линейные полиномы

u= +

v= + (10.4)

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

Система уравнений для перемещений u узлов i, j, k:

= + ;

= + ; (10.5)

= + ;

позволяет определить коэффициенты , , по правилу Крамера через определители:

; ; ,

где ∆ — площадь треугольника ijk, удвоенное значение которой равно определителю

2∆=det = . (10.6)

После подстановки , , в первое уравнение (10.4), получим

u= , (10.7)

где коэффициенты определяются через координаты узлов КЭ:

; ; = . (10.8)

Остальные коэффициенты получаются циклической перестановкой индексов i, j, k.

Аналогично можно представить перемещение v, пользуясь вторым уравнением (10.4):

v= . (10.9)

Вектор перемещений внутри КЭ «e»принимает вид

= = , (10.10)

где функции формы , , , входящие в матрицу формы , определяются зависимостями:

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

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

= = . (10.11)

Подставляя (10.10) в (10.11) получим

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

Матрица напряжений.Для упругого изотропного материала физические соотношения т.е. зависимости между напряжениями и деформациями, в данном случае закон Гука, линейны и определяются уравнениями:

где ̶ матрица упругости, содержащая упругие постоянные материала. Для изотропного материала в случае обобщенного плоского напряженного состояния ОПНС :

где 𝐸 ‒модуль упругости; 𝜇 ‒ коэффициент Пуассона.

Используя (10.12) и (10.13), напряжения в КЭ 𝑒можно выразить через перемещения его узлов

где матрица называется матрицей напряжений:

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

Матрица жесткости конечного элемента 𝑒 (𝑖𝑗𝑘)устанавливает соответствие между узловыми силами и перемещениями по направлению этих сил. Матрица осредненно характеризует жесткостные свойства среды в объеме -го конечного элемента, а вектор ‒ приходящееся на этот элемент внешнее воздействие. Известно, что решение плоской задачи эквивалентно отысканию функций и , удовлетворяющих на контуре краевым условиям и минимизирующих функционал

где энергия деформации 𝑈 определяется выражением (в соответствии с формулой Клапейрона):

Энергия деформации 𝑈 для КЭ 𝑒, с учетом матричного представления и полученных ранее соотношений

Напомним, что транспонирование произведения двух матриц выполняется по формуле (лк.1 ЧМ ч.1):

Работа 𝐴, совершаемая приложенными силами, может быть разделена на три части:

‒ работа совершаемая внешними сосредоточенными силами в узлах;

‒работа , совершаемая распределенными поверхностными силами приложенными на гранях КЭ;

‒работа , выполняемая объемными силами 𝑋, 𝑌.

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

Поверхностные и объемные силы можно также привести к узловым. Тогда

𝐴 = = .

При этом выражение для полной потенциальной энергии для КЭ принимает вид:

Чтобы согласно вариационному принципу Лагранжа минимизировать величину , необходимо продифференцировать выражение (10.18) по и результат приравнять нулю:

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

Теперь матрица жесткости элемента 𝑖𝑗𝑘 определяется так:

Отметим единообразие выражений для матриц жесткости стержневых и плоских (10.20) КЭ.

Если толщина элемента 𝑒 постоянна, то:

= h∆ ,(10.21)

где: Δ ‒ площадь КЭ; ‒ его объем.

Матрица жесткости является симметричной и имеет блочную структуру:

=

где ‒ квадратная симметричная подматрица порядка 𝑚 (𝑚 ‒ число компонент перемещений в узле), элементами которой являются реактивные усилия в связях, наложенных на узел 𝑟, вызванные единичными перемещениями и узла 𝑠.

Реактивные усилия в связях, наложенных на узлы 𝑖, 𝑗, 𝑘, вызванные перемещениями всех узлов, равны:

так, например, — реактивные усилия в связях, наложенных на узел 𝑖, вызванные перемещениями всех узлов

При ОПНС подматрицы имеют второй порядок и определяются следующим образом:

где и ‒ блоки градиентной матрицы :

;

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

Узловые силы, связанные с объемными силами X и Y (проекции их интенсивностей на оси x и y), определяются выражениями:

dxdy ,

= — , и т.д.

Т.к. функции формы матрицы зависят от переменных x и y, надо интегрировать выражение . Если за начало координат выбрать центр тяжести КЭ, то

= — ∆/3 (10.26)

Узловые силы, вызванные поверхностными нагрузками и ,равны:

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

10.3. Построение глобальной матрицы жесткости. Поскольку полная потенциальная энергия (функционал) для рассматриваемой области может быть представлен суммой всех КЭ

то, минимизируя этот функционал, получим следующую общую систему уравнений:

,

где – вектор, содержащий перемещения всех узлов рассматриваемой области :

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

Пусть – типовой узел расчетной схемы области .

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

Для показанного на рис. 1 конечного элемента , . Вектор реакций можно записать в виде:

Из условий равновесия узла следует:

где – вектор внешних сил, приложенных к узлу , который может формироватся за счет объемных сил, теплового воздействия , а для контурных КЭ – и поверхностных сил.

Вектор можно представить и так:

где – число узлов конечных элементов ;

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

Из подстановки (10.32) в (10.31) следует

где – матрица жесткости всего диска.

– вектор узловых сил, эквивалентных внешней нагрузке.

Линейные алгебраические уравнения типа (10.33) составляются для всех узлов расчетной схемы области и образуют систему уравнений для определения узловых перемещений .

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

Реакция в наложенной связи, исключавшей перемещение можно найти по формуле

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

Метод конечных элементов. Метод конечных разностей

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

Рассмотрим основные понятия метода конечных разностей (МКР) на примере наиболее простого уравнения – уравнения переноса величины 𝑢(𝑥,𝑡):

Физический смысл величины 𝑢 может быть различным. Например, это может быть давление в звуковой волне или электрический сигнал.

Физический смысл уравнения (1) заключается в следующем. В начальный момент времени (𝑡=0) функция 𝑢(𝑥,0) имеет некоторый профиль 𝑔(𝑥), т.е. это начальное условие задачи:

Точное решение уравнения дается выражением 𝑢(𝑥,𝑡)=𝑔(𝑥−𝑐𝑡), т.е. решение представляет собой ту же начальную функцию 𝑔(𝑥), но смещенную относительно начального положения на расстояние 𝑐𝑡:

Таким образом, уравнение (1) описывает перемещение начального профиля 𝑔(𝑥) вправо по направлению оси 𝑥 со скоростью 𝑐. Поэтому уравнение называется уравнением переноса.

Применение МКР для решения уравнения переноса

Первым шагом для применения МКР является введение разностной сетки, показанной на рисунке:

Сетка может быть равномерной (шаги сетки одинаковые) и неравномерной (шаги сетки различны).

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

Индекс 𝑘 относится ко времени, индекс 𝑛 – к координате.

Это называется сеточной функцией. Очевидно, что эта функция дискретна, т.к. область ее значений – точки в узлах сетки.

Производные на сетке можно аппроксимировать так (№1):

где h – шаг по координате.

Аналогично поступают с производными по времени, где τ – шаг по времени.

Тогда разностную схему для уравнения переноса можно записать следующим образом:

Из этой разностной схемы легко выразить значение 𝑢 на 𝑘+1 слое по времени:

И таким образом вычислить решение итеративно. Поэтому такая схема называется явной.

С другой стороны, можно аппроксимироватьпроизводные иначеи получить другую разностную схему (№2):

Она содержит значения 𝑢 на 𝑘+1 слое в трех местах, поэтому отсюда его явно выразить нельзя. В результате приходится записывать сразу все уравнения на всей разностной сетке и решать систему уравнений на каждом шаге по времени. Здесь используются методы, рассмотренные нами ранее для решения уравнений и систем уравнений.

Наконец, можно рассмотреть и такую схему (№3):

Откуда можно выразить 𝑢_𝑛 𝑘 +1 .

Метод конечных объемов

Метод конечных объемов подобен методу конечных разностей

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

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

Метод конечных элементов (МКЭ, Finite Elements Method, FEM)– это численный метод решения дифференциальных уравнений с частными производными, а также интегральных уравнений, возникающих при решении задач прикладной физики. Метод широко используется для решения задач механики деформируемого твёрдого тела, теплообмена, гидродинамики и электродинамики.

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

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

Рассмотрим применение МКЭ чуть подробней. В общем случае исследуемая непрерывная величина заранее неизвестна, и нужно определить значения этой величины в некоторых внутренних точках области. Дискретную модель, однако, очень легко построить, если сначала предположить, что числовые значения этой величины в каждой внутренней точке области известны. После этого можно перейти к общему случаю. Итак, при построении дискретной модели непрерывной величины поступают следующим образом:

1. В рассматриваемой области фиксируется конечное число точек. Эти точки называются узловыми точками или просто узлами.

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

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

4. Непрерывная величина аппроксимируется на каждом элементе полиномом, который определяется с помощью узловых значений этой величины. Для каждого элемента определяется свой полином, но полиномы подбираются таким образом, чтобы сохранялась непрерывность величины вдоль границ элемента. Основная идея МКЭ наглядно показана на одномерном примере заданного распределения температуры в стержне рис. 1

Рассматривается непрерывная величина T(x), область определения которой – отрезок OL вдоль оси х. Фиксированы и пронумерованы пять точек на оси x (рис. 2, а). Это узловые точки. Совсем не обязательно располагать узловые точки на равном расстоянии друг от друга. Очевидно, можно ввести в рассмотрение и более пяти точек, но этих пяти вполне достаточно, чтобы проиллюстрировать основную идею метода. Значения T(x) в данном случае известны в каждой узловой точке. Эти фиксированные значения представлены графически на рис. 2, б и обозначены в соответствии с номерами узловых точек через T₁, T₂, T₃, T₄, T₅.

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

Окончательная аппроксимация Т(х) будет состоять из четырех кусочно-линейных функций, каждая из которых определена на отдельном элементе (рис. 4).

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

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