Численные методы решения задач тепломассообмена и гидродинамики

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

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

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

8.2 Метод конечных разностей: основные понятия

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

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

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

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

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

- каким образом выбрать сетку;

- как построить разностную схему;

- определить точность, с которой разностная схема аппроксимирует исходную задачу;

- проверить устойчивость разностной схемы;

- выяснить скорость сходимости решения разностной задачи к решению исходной.

8.2.1 Построение сеток. Сеточные функции и сеточные аналоги норм

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

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

Примеры сеток.

А) В простейшем случае одномерной задачи отрезок [0, l] разобьем на N равных частей точками хk=kh, k=0,1,…, N

Расстояние между узлами h=(хk+1 - хk) называется шагом сетки. В рассматриваемом случае h = l/N = const, поэтому множество узлов хk, k=0,1,…, N представляет собой равномерную сетку на отрезке 0 ≤ х ≤ l и

 
  Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

обозначается Ωk = { хk=kh, k=0,1,…, N }.

0 h хk = kh хk+1/2 = хk +h /2 l

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

0 1 2 … k-1 k k+1…. N-1 N x

Рисунок 8.1 – Одномерная равномерная сетка на отрезке

В качестве области определения сеточных функций кроме узлов, называемых еще целыми точками, часто используют полуцелые точки хk+1/2 = хk +h /2 (на рисунке отмечены крестиком).

Если отрезок [0, l] разобьем на N частей произвольно взятыми точками 0 <х123…< хk< хk+1< …< хN-1< l, причем х0=0, хN=l, то получим неравномерную сетку с шагом hk= хk- хk-1, зависящем от номера узла k.

Б) Сетка на плоскости.

Пусть Ω = {0 ≤ х ≤ d, 0 ≤ y ≤ b } – прямоугольник. Отрезки [0, d] и [0, b] разобьем соответственно на N и M частей и через точки хk=kh, k=0,1,…, N, h=d/N, yj=jρ, j=0,1,…, M, ρ=b/M, проведем прямые, параллельные координатным осям. Множество точек (xk,yj) образует сетку в прямоугольнике Ω. Полученная сетка равномерна по каждой переменной х и у. Если h #ρ – сетка прямоугольная, при h = ρ – квадратная.

Если построить сетку неравномерной хотя бы по одной координате, то сетка будет неравномерной.

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru
Рисунок 8.2 - Равномерная сетка на плоскости

Для решения одномерной по пространственным координатам нестационарной задачи используют произведение сеток

           
  Численные методы решения задач тепломассообмена и гидродинамики - student2.ru   Численные методы решения задач тепломассообмена и гидродинамики - student2.ru   Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

Ωh х Ωη={(хk, τi ), хk+1= хk + hk, τi+1 = τi + ηi , k=0,1,…, N, i=0,1, …, M, х0=0, хN=l, τ0=0, τM=t}, представляющих собой пространственно-временную разностную сетку.

Совокупность узлов сетки, лежащих на τ = τi, называют i-м слоем.

Для простых областей всегда можно ввести такую сетку, чтобы «крайние» естественнее узлы сетки попадали на границу области. Эти узлы называются граничными, а остальные – внутренними. Граничные условия задачи следует задавать именно в этих граничных узлах.

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

- ввести дополнительные узлы в точках пересечения линий сетки с границей и в них задавать граничные условия;

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

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

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

Для определения степени близости этих функций обычно поступают так. Осуществляют переход от непрерывных функций к сеточным по правилу: значение сеточной функции в узле равно значению непрерывной функции в этой же точке. Например, в узле (хk, τi) сетки одномерной нестационарной задачи Tki=T(хk, τi) пространственнее узлы сетки обозначают подстрочными индексами, временные – надстрочными. В этом случае говорят так, что сеточная функция Tki является проекций функции Т на пространство сеточных функций, либо Тh. Для оценки близости функций u и Т рассматривается величина ||u -Тh||, где ||·|| некоторая норма в пространстве сеточных функций. Чаще всего ||u||С = mах ||u|| для хk є Ωk – сеточный аналог чебышевской нормы в пространстве непрерывных функций [].

8.2.2 Построение разностных схем. Порядок аппроксимации

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

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

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

- метод формальной замены производных конечно-разностными выражениями;

- метод интегральных тождеств (интегр-интерполяционный метод);

- вариационные методы построения разностных схем;

- метод неопределенных коэффициентов.

8.2.3 Метод формальной замены производных конечно-разностными выражениями

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

Рассмотрим возможные способы аппроксимации дифференциального оператора вида

А[T] =dT/dx, (8.1)

определенного на множестве непрерывных в области Ω={d<x<b} функций, имеющих ограниченные производные до третьего порядка включительно. Пусть Ωh = {хk=kh, 0<k<N, h=(b-d)/N} – равномерная сетка на отрезке Ω. Наиболее естественным способом замены производной будет

dT/dx = lim [ T(x+h) - T(x)]/h

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru h 0

Если зафиксировать h в этом равенстве, то получим приближенную формулу для первой производной через конечные разности

dT/dx ≈ [ T(x+h) - T(x)]/h

или в k-м узле имеем правое разностное отношение

Tx,k = (Tk+1 -Tk)/h . (8.2)

Аналогично вводится левое разностное отношение

T*x,k = (Tk -Tk-1)/h . (8.3)

Можно рассмотреть линейную комбинацию (8.2) и (8.3):

σ Tx,k+ (1 - σ) T*x,k , (8.4)

где σ – любое вещественное число.

При σ=1/2 получим центральное разностное отношение:

Tx,k 0 = (Tk+1 -Tk-1)/2h. (8.5)

При замене оператора А [T] разностными выражениями (8.2) –(8.5) допускается погрешность аппроксимации 0(hn), где n=1 для (8.2) – (8.4) и

n =2 для (8.5).

Таким образом, выбор шаблона существенно влияет на свойства разностного оператора.

Для аппроксимации второй производной

А[T] =d2T/dx2

выбираем трехточечный шаблон (хk-1, xk, xk+1).

Тогда в k-м узле получим разностный оператор

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru Txx,k= (Tk+1 -2Tk + Tk-1)/h2 . (8.6)

Порядок аппроксимации в этом случае равен 2, погрешность аппроксимации 0(h2).

Рассмотрим оператор одномерного уравнения теплопроводности

А[T] = dT/dτ - аd2T/dx2 (8.7)

в области Ω = {(0<х<d, 0<τ<t)}, предполагая, что функция температуры непрерывна со своими производными до четвертого порядка по переменной х и до второго по - τ.

Область непрерывного изменения аргументов заменим сеточной областью Ω={(хk, τi ), хk =k·h, τi =i·η , 0 <k< N, 1< i < M, η=t/M }.

Аппроксимируем производную по времени правым разностным соотношением Tiτ = (Tk i+1 - Tk I)/ η, а для второй производной по переменной х можно записать отношение (8.6) на временном слое i:

 
  Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

Txx,k i = (Tk+1i -2Tk i + Tk-1 i)/h2 ,

или на временном слое i+1:

 
  Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

Txx,k i+1 = (Tk+1i+1 -2Tk i+1 + Tk-1 i+1)/h2 .

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru В соответствии с этим можно рассмотреть две различные аппроксимации оператора (8.7):

А[T] = Tiτ - а Txx,k i (8.8)

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru А[T] = Tiτ - а Txx,k i+1 (8.9)

на шаблонах явной и неявной схем (рис 8.3).

    Численные методы решения задач тепломассообмена и гидродинамики - student2.ru     Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

а) б)

Рисунок 8.3 - Шаблоны явной (а ) и неявной (б) схем

Погрешность локальной аппроксимации оператора (8..7) разностными операторами (8.8) и (8.9) равна 0(η + h2).

Оператор (8.7) можно аппроксимировать и на шеститочечном шаблоне Кранка-Николсона

Численные методы решения задач тепломассообмена и гидродинамики - student2.ru

Рисунок 8.4 – Шаблон неявной схемы Кранка-Николсона

Погрешность локальной аппроксимации оператора (8.7) разностным оператором при этом равна 0(η2 + h2).

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

Компьютерная реализация методов описана в [9-10], [23-26].

Литература: [1-4],[6], [9-10], [14-17], [23-24].

Лекция 9

Наши рекомендации