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

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

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

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

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 и

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

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

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

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

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

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

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

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

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

Рисунок 8.2 - Равномерная сетка на плоскости

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

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

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

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

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

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

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

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

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

Для определения степени близости этих функций обычно поступают так. Осуществляют переход от непрерывных функций к сеточным по правилу: значение сеточной функции в узле равно значению непрерывной функции в этой же точке. Например, в узле (х_k, τ_i) сетки одномерной нестационарной задачи T_kⁱ=T(х_k, τ_i) пространственнее узлы сетки обозначают подстрочными индексами, временные – надстрочными. В этом случае говорят так, что сеточная функция T_kⁱ является проекций функции Т на пространство сеточных функций, либо Т_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

h 0

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

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

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

T_x,k = (T_k+1 -T_k)/h . (8.2)

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

T*_x,k = (T_k -T_k-1)/h . (8.3)

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

σ T_x,k+ (1 - σ) T*_x,k , (8.4)

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

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

T_x,k ⁰ = (T_k+1 -T_k-1)/2h. (8.5)

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

n =2 для (8.5).

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

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

А[T] =d²T/dx²

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

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

T_xx,k= (T_k+1 -2T_k + T_k-1)/h² . (8.6)

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

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

А[T] = dT/dτ - аd²T/dx² (8.7)

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

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

Аппроксимируем производную по времени правым разностным соотношением Tⁱ_τ = (T_k ⁱ⁺¹ - T_k ^I)/ η, а для второй производной по переменной х можно записать отношение (8.6) на временном слое i:

T_xx,k ⁱ = (T_k+1ⁱ -2T_k ⁱ + T_k-1 ⁱ)/h² ,

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

T_xx,k ⁱ⁺¹ = (T_k+1ⁱ⁺¹ -2T_k ⁱ⁺¹ + T_k-1 ⁱ⁺¹)/h² .

В соответствии с этим можно рассмотреть две различные аппроксимации оператора (8.7):

А[T]_hη = Tⁱ_τ - а T_xx,k ⁱ (8.8)

А[T]_hη = Tⁱ_τ - а T_xx,k ⁱ⁺¹ (8.9)

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

а) б)

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

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

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

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

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

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

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

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

Лекция 9

Тема: Численные методы решения задач оптимизации

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