Конечно разностный аналог уравнения теплопроводности

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

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

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

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

1 Теоретическая часть

1.1 Постановка задач для уравнений параболического типа

Классическим примером уравнения параболического типа является уравнение теплопроводности (диффузии). В одномерном по пространству случае однородное (без источников энергии) уравнение теплопроводности имеет вид

(1)

Если на границах и заданы значения искомой функции в виде

, , (2)

, , (3)

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

, , (4)

то задачу (1)-(4) называют первой начально-краевой задачей для уравнения теплопроводности (1).

В терминах теории теплообмена — распределение температуры в пространственно-временной области

a 2 — коэффициент температуропроводности, а (2), (3) с помощью функций , задают температуру на границах и .

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

, , (5)

, , (6)

т.е. граничные условия второго рода, то задачу (1), (5), (6), (4) называют второй начально-краевой задачей для уравнения теплопроводности (1). В терминах теории теплообмена на границах в этом случае заданы тепловые потоки.

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

, , (7)

, , (8)

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

Для пространственных задач теплопроводности в области первая начально-краевая задача имеет вид

(9)

Аналогично ставится вторая и третья начально-краевые задачи для пространственного уравнения (9). На практике часто ставятся начально-краевые задачи теплопроводности со смешанными краевыми условиями, когда на границах задаются граничные условия различных родов.

1.2 Основные определения и конечно-разностные схемы

Основные определения, связанные с методом конечных разностей, рассмотрим на примере конечно-разностного решения первой начально-краевой задачи для уравнения теплопроводности (1)-(4).

Согласно методу сеток в плоской области D строится сеточная область D_h , состоящая из одинаковых ячеек. При этом область D_h должна как можно лучше приближать область D . Сеточная область (то есть сетка) D_h состоит из изолированных точек, которые называются узлами сетки. Число узлов будет характеризоваться основными размерами сетки h : чем меньше h , тем больше узлов содержит сетка. Узел сетки называется внутренним, если он принадлежит области D , а все соседние узлы принадлежат сетке D_h . В противном случае он называется граничным. Совокупность граничных узлов образует границу сеточной области Г_h .

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

Нанесем на пространственно-временную область , конечно разностную сетку ω_h,τ :

(10)

с пространственным шагом h = l / N и шагом по времени τ=T/K.

Рисунок 1 – Конечно-разностная сетка

Введем два временных слоя : нижний ,на котором распределение искомой функции u ( x_j , t k ) , , известно (при к = 0 распределение определяется начальным условием (4)u ( x_j , t k )=ψ( x_j ) ), и верхний временной слой t k +1 =( k +1) τ , на котором распределение искомой функции u ( x_j , t k +1 ) , .

Сеточной функцией задачи (1)-(4) называют однозначное отображение целых аргументов j , k в значения функции .

На введенной сетке вводят сеточные функции , первая из которых известна, вторая подлежит определению. Для определения в задаче (1)-(4) заменяют (аппроксимируют) дифференциальные операторы отношением конечных разностей (более подробно это рассматривают в разделах численных методов «Численное дифференцирование»), получают

, (11)

, (12)

Подставляя (11), (12) в задачу (1)-(4), получим явную конечно-разностную схему для этой задачи в форме

(13)

В каждом уравнении этой задачи все значения сеточной функции известны, за исключением одного, , которое может быть определено явно из соотношений (13). В соотношения (13) краевые условия входят при значениях j =1 и j = N — l , a начальное условие – при k = 0.

Если в (12) дифференциальный оператор по пространственной переменной аппроксимировать отношением конечных разностей на верхнем временном слое:

, (14)

то после подстановки (11), (14) в задачу (1)-(4) получим неявную конечно-разностную схему для этой задачи:

(15)

Теперь сеточную функцию на верхнем временном слое можно получить из решения (15) с трехдиагональной матрицей. Эта СЛАУ в форме, пригодной для использования метода прогонки, имеет вид

;

;

, ;

;

;

;

.

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

Рисунок 2 — Шаблон явной конечно-разностной схемы для уравнения теплопроводности

Рисунок 3 — Шаблон неявной конечно-разностной схемы для уравнения теплопроводности

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

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

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

Вопрос устойчивости будет рассмотрен далее.

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

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

1.4 Устойчивость. Исследование устойчивости методом гармонического анализа

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

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

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

, (16)

где λ _n – собственные значения

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

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

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

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

(17)

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

(18)

где – размерный и постоянный сомножитель амплитудной части,

k – показатель степени (соответствующий номеру временного слоя) сомножителя, зависящего от времени.

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

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

(19)

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

, (20)

т. е. условие (20) является необходимым условием устойчивости.

1.5 Схема Кранка-Николсона

параболическое дифференциальное уравнение конечная разность

Явная конечно разностная схема, записанная в форме

(21)

обладает тем достоинством, что решение на верхнем временном слое t_k+l получается сразу (без решения СЛАУ) по значениям сеточной функции на нижнем временном слое t k , где решение известно (при k = 0 значения сеточной функции формируются из начального условия). Но эта же схема обладает существенным недостатком, поскольку она является условно устойчивой. С другой стороны, неявная конечно-разностная схема, записанная форме

(22)

приводит к необходимости решать СЛАУ, но зато эта схема абсолютно устойчива.

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

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

На убывающем решении картина изменяется противоположным образом: явная конечно-разностная схема завышает решения, а неявная — занижает (Рисунок 4).

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

Рисунок 4 – Двусторонний метод аппроксимации

Проведенный анализ дал блестящий пример так называемых двусторонних методов, исследованных В. К. Саульевым

Рассмотрим неявно-явную схему с весами для простейшего уравнения теплопроводности:

(23)

где θ – вес неявной части конечно-разностной схемы,

θ -1 – вес для явной части

Причем . При θ=1 имеем полностью неявную схему, при θ=0 – полностью явную схему, а при θ=1/2 – схему Кранка-Николсона .

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

,

(24)

причем правое неравенство выполнено всегда.

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

(25)

являющаяся условием устойчивости неявно-явной схемы с весами (23), когда вес находится в пределах .

Таким образом, неявно-явная схема с весами абсолютно устойчива при и условно устойчива с условием (25) при .

Рассмотрим порядок аппроксимации неявно-явной схемы с весами, для чего разложим в ряд Тейлора в окрестности узла (x j ,t_k ) на точном решении значения сеточных функций по переменной t , , по переменной х и полученные разложения подставим в (23):

В этом выражении дифференциальный оператор от квадратной скобки в соответствии с дифференциальным уравнением равен дифференциальному оператору , в соответствии с чем вышеприведенное равенство приобретает вид

После упрощения получаем

,

откуда видно, что для схемы Кранка-Николсона (θ = 1/2) порядок аппроксимации схемы (23) составляет , т.е. на один порядок по времени выше, чем для обычных явных или неявных схем. Таким образом, схема Кранка-Николсона при θ = 1/2 абсолютно устойчива и имеет второй порядок аппроксимации по времени и пространственной переменной х .

Используем в уравнение (23) подстановку r= a 2 k / h 2 . Но в то же время его нужно решить для трех «еще не вычисленных» значений , , и . Это возможно, если все значения перенести в левую часть уравнения. Затем упорядочим члены уравнения (23) и в результате получим неявную разностную формулу

(26)

для i=2,3,…, n-1 . Члены в правой части формулы (26) известны. Таким образом, формула (26) имеет вид линейной трехдиагональной системы АХ=В. Шесть точек, используемых в формуле Кранка-Николсона (26), вместе с промежуточной точкой решетки, на которой основаны численные приближения, показаны на рисунке 5.

Рисунок 5 – Шаблон (схема) метода Кранка-Николсона

Иногда в формуле (26) используется значение r=1 . В этом случае приращение по оси t равно , формула (26) упрощается и принимает вид

, (27)

для i=2,3,…, n-1 . Граничные условия используются в первом и последнем уравнениях (т. е. в и соответственно).

Уравнения (27) особенно привлекательны при записи в форме трехдиагональной матрицы АХ = В.

Если метод Кранка-Николсона реализуется на компьютере, то линейную систему АХ = В можно решить либо прямым методом, либо итерационным.

2. Практическая часть

2.1 Постановка задачи

Используем метод Кранка-Николсона, чтобы решить уравнение

,

с начальным условием

,

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

2.2 Решение в ППП MatLab

Решение будем искать в ППП MatLab 7. Создадим четыре выполняемых m-фала: crnich.m – файл-функция с реализацией метода Кранка-Николсона; trisys.m – файл-функция метода прогонки; f.m – файл-функция задающая начальное условие задачи; fе.m – файл-функция задающая функцию определяющую точное решение задачи(найдена аналитическим путем). Листинги программ представлены в приложении А.

Для простоты возьмем шаг Δх = h = 0,1 и Δ t = к = 0,01 . Таким образом, соотношение r =1. Пусть решетка имеет n=11 столбцов в ширину и m=11 рядов в высоту.

2.3 Анализ результатов

Решения для данных параметров отразим в таблице 1. Трехмерное изображение данных из таблицы покажем на рисунке 5.

Таблица 1 – Значения u(х _i , t_i ), полученные методом Кранка-Николсона

Величины, полученные методом Кранка-Николсона, достаточно близки к

аналитическому решению u(x,t) = sin(πx)e -π2 t + sin(3πx)e -9π2 t , истинные значения для последнего представлены в таблице 2

Максимальная погрешность для данных параметров равна 0,005

Таблица 2 – точные значения u(х _i , t_i ), при t=0.1

Рисунок 5 –Решениеu= u(х _i , t_i ), для метода Кранка-Николсона

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

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

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

1 Березин И.С., Жидков Н.П. Методы вычислений. Т.2. – М.: Физматгиз, 1962.

2 Мэтьюз, Джон, Г., Финк, Куртис, Д. Численные методы. Использование MATLAB, 3-е издание.— М. : Вильяме, 2001. — 720 с

3 Тихонов А.Н., Самарский А.А. Уравнения математической физики. – М.: Наука, 1972.

4 Формалев В.Ф., Ревизников Д.Л. Численные методы. – М.: ФИЗМАТЛИТ, 2004. — 400 с.

5 Пирумов У.Г. Численные методы. – М.: Издательство МАИ, 1998.

6 Калиткин Н.Н. Численные методы. – М.: Наука, 1976.

Листинг программы для расчета по методу Кранка-Николсона

Конечно-разностные аппроксимации производных (стр. 1 )

Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6 7

Конечно-разностные аппроксимации производных

Конечно-разностные аппроксимации производных (конечные разности) — способ приближенного вычисления частных производных

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

Или более коротко с использованием индексов точек:

Отсюда , где — остаток.

Отбрасывая остаток можно получить правую разность:

Погрешность такой аппроксимации определяется старшим членом в отброшенном остатке и в данном случае этот член содержит в первой степени.

Аналогичным образом, разлагая в ряд функцию можно получить:

Получим новую аппроксимацию первой производной:

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

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

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

Отброшенный остаток будет содержать член с во второй степени (после деления на )

Исходя из определения, можно получить выражения для третьей, четвертой и более старших разностей:

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

правая разность по оси : ; правая разность по оси : ; левая разность по оси : ; левая разность по оси : ; центральная разность по оси : ; центральная разность по оси : ; вторая разность по оси : ; вторая разность по оси : .

Смешанная производная может быть получена следующим образом:

Алгоритм решения стационарных краевых задач методом конечных разностей

Метод конечных разностей — универсальный сеточный численный метод решения задач микроуровня.

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

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

Решение линейных одномерных стационарных краевых задач с помощью МКР

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

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

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

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

для узла 1:
для узла 2:

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

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

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

Пусть и .

Запишем разностные аналоги для внутренних узлов сетки:

для узла 1:
для узла 2

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

Проще всего воспользоваться левой разностью:

Решая эту систему уравнений, получим , , .

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

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

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

для узла 1:
для узла 2:
для узла 3:
граничное условие второго рода:

За повышение точности пришлось заплатить увеличением размерности системы конечно-разностных уравнений.

Решение нелинейных одномерных стационарных краевых задач с помощью МКР

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

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

где — коэффициент теплопроводности.

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

Предположим, что толщина стекла . Температура в комнате , на улице —

Тепловой поток на улицу пропорционален градиенту температуры, то есть .

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

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

для узла 1:
для узла 2:
для узла 3:

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

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

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

для узла 1:
для узла 2:

Данную систему придется решать итерационными методами.

Решение нестационарных одномерных задач с помощью МКР

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

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

Пусть при записи разностей нижний индекс соответствует оси , а верхний — оси времени.

Первый вариант разностного уравнения, апроксимирующего исходное (1):

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

Второй вариант разностного уравнения, апроксимирующего исходное (1):

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

Графическое изображение разностных уравнений получило название шаблонов решения сответствующих задач. В данном случае на рис. 1,а представлен шаблон явной разностной схемы, а на рис. 1,б — неявной.

Рис. 1. Шаблоны явной и неявной разностной схемы

Использование шеститочечного шаблона применено в схеме Кранка-Николсона:

В общем случае использования шеститочечного шаблона, имеем схему с весами:

которая при является неявной.

Примеры решения нестационарных задач с помощью МКР

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

Решим задачу с помощью явной разностной схемы.

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

Пусть , выберем значения шага по оси и значение шага по оси времени .

Наносим на объект равномерную сетку по оси , как это показано на рис. 1.

Записываем явную разностную схему для узла 1:

где — граничное условие, — начальные условия, отсюда .

Записываем явную разностную схему для узла 2:

где — граничное условие, — начальные условия, отсюда .

Таким образом найдено температурное поле в момент времени .

Аналогично для момента времени :

Для момента времени :

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

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

Как нетрудно проверить, условие не было выполнено. Чтобы удостовериться в работоспособности явной разностной схемы, повторим вычисления для :

для момента времени ; для момента времени для момента времени .

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

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

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

Рассмотрим решение задачи явной разностной схемой с граничными условиями второго рода (типа Неймана).

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

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

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

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

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

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

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

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

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

При переходе к тепловой цепочке уравнение проводимости между

ее узлами может быть записано в виде

где R_n-1,_n— сосредоточенное термическое сопротивление между узлами n- 1 и n; t_n-1 и t_n-температуры в узлах тепловой цепочки, где сосредоточены теплоемкости.

Уравнение для тепловой цепочки справедливо как для стационарных, так и нестационарных условий.

Закон аккумуляции тепла устанавливает, что при ращение количества тепла dQ, аккумулированного слоем dx, пропорционально приращению во времени его температуры

где ср – объемная теплоемкость материала.

Изменение количества аккумулированного тепла ΔQ для элементарного слоя толщиной Δх при изменении во времени z его средней температуры на Δ_zt равно

Для тепловой цепочки уравнение аккумуляции тепла может быть

записано в виде

где С = срΔх — сосредоточенная тепловая емкость элементарного слоя; Δ_zt -изменение во времени (z) температуры в центре элементарного слоя в сечении расположения сосредоточенной емкости.

Составим уравнение теплового баланса элементарного слоя n при распределении температур в сечении, отмеченном на рисунке tox. Слой n обменивается теплом с соседними элементарными слоями и согласно закону проводимости за время Δz он получит от слоя n + 1 количество тепла

и отдаст слою n -1 количество тепла

Разность ΔQn между количествами тепла, определенными этими уравнениями, будет аккумулирована слоем n и повысит его среднюю температуру на Δ_ztn .

Уравнение теплового баланса слоя n можно написать в виде

которое после преобразований может быть записано

где

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

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

Применительно к тепловой цепочке уравнения теплопроводности в конечных

Множитель в виде комплекса величин в левой части этого уравнения является обратной величиной критерия гомохронности (Фурье) процесса, написанного для элементарного слоя Δх и расчетного интервала времени Δz. После подстановки значений этот множитель можно преобразовать и заменить обозначением критерия Фурье:

а — коэффициент температуропроводности.

Тогда уравнение теплопроводности в конечных разностях принимает вид

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

Дата добавления: 2014-12-03 ; просмотров: 12 ; Нарушение авторских прав