МОДЕЛИРОВАНИЕ ПРОЦЕССА ТЕПЛОПРОВОДНОСТИ
То, что тела могут проводить тепло, общеизвестно. Если один из концов длинного стержня поместить в костер, то, если стержень сделан не из горючего или легко плавящегося материала, другой конец через некоторое время тоже нагреется; как быстро и насколько - зависит от материала, размеров стержня и других факторов. Процесс теплопроводности - один из, так называемых, процессов тепломассопереноса, играющих огромную роль в природе и в технике. Другие процессы такого рода - диффузия, благодаря которой смешиваются разные жидкости или газы, процессы гидро- и аэродинамики (т.е. переноса (движения) жидкостей и газов).
Хотя каждый из таких процессов имеет собственные закономерности, между ними много общего. Эти процессы происходят в сплошной среде, о которой шла речь выше; при их математическом моделировании используется один и тот же математический аппарат-дифференциальные уравнения в частных производных.
Ограничимся одной из самых простых задач данного класса - переносом тепла в однородном стержне. Рассмотрим линейный стержень, боковая поверхность которого не проводит тепла (теплонзолирована). Если в начальный момент стержень неравномерно нагрет, то в нем будет происходить перераспределение тепла; при отсутствии внутренних источников тепла его температура, в конце концов, выровняется.
Поскольку стержень линеен и однороден, то распределение температуры в пространстве характеризуется одной координатой x.
Температура (обозначим ее u) зависит от х; кроме того, она может меняться со временем, т.е. является функций двух переменных и(х, t). Изменение этой функции вдоль стержня, «скорость» которого определяется производной пол x, и изменение ее со временем, скорость которого определяется производной по t, взаимосвязаны и, как будет показано ниже, входят в одно уравнение.
Уравнение теплопроводности.
Получим уравнение, описывающее процесс изменения температуры в стержне. Фиксируем некоторую точку x0
(рис. 7.29) и выделим около нее малый участок стержня длиной ?x.
Искомое уравнение есть по существу уравнение теплового баланса (т.е. сохранения энергии): изменение количества тепла в избранном участке стержня за счет притока и (или) оттока его через два сечения приведет к нагреванию или охлаждению этого участка в соответствии с его теплоемкостью. Выразим все это математическим языком.
Рис. 7.29. Участок линейного стержня
Количество тепла, проходящее через поперечное сечение стержня в точке x0
за время ?t, пропорционально площади поперечного сечения S,
градиенту температуры
Вводя коэффициент пропорциональности k, называемый коэффициентом теплопроводности, получаем
Значение k определяется материалом стержня и для нескольких материалов приведено в табл. 7.6 (в единицах системы СИ:
Таким образом, различия в теплопроводности разных материалов огромны.
Рис. 7.30. Поток тепла через участок стержня длиной ?х
Теперь запишем количество тепла, проходящее через сечение в точке х = x0 + ?x:. Оно определяется, естественно, той же формулой:
с условием, что производная
Таблица 7.6
Значение коэффициента теплопроводности для некоторых материалов
|
Медь |
384 |
Лед (0° С) |
2,23 |
Асбест |
0,4 - 0,8 |
|
Алюминий |
209 |
Бетон |
0,7 - 0,2 |
Дерево |
0,1 - 0,2 |
|
Сталь |
47 |
Кирпич |
0,7 |
Воздух |
0,034 |
?x,
в силу чего
Если через сечения х = х0
и х = x0 + ?x за время ?t прошло разное количество тепла, то та его часть, которая пошла на нагревание (или, в зависимости от знака, на охлаждение) этого участка стержня, есть
Пусть за то же время температура участка изменилась на ?u; как известно, это связано с изменением ?Q соотношением ?Q
= mc?u, где т - масса, с - удельная теплоемкость. Приравняем два выражения для ?Q:
Поскольку массу можно представить как т = ?•S•?x (? - плотность вещества), то, поделив обе части уравнения на ?t и перейдя к пределу при ?t > 0, получим
Это - основное уравнение теплопроводности для однородного стержня. Как следует из процедуры вывода, это уравнение локально, т.е. в данный момент времени и в данной точке выражает закон сохранения энергии.
В уравнение (7.48) входят три постоянные, характеризующие вещество. Удобно объединить их в одну, переписав уравнение в виде
где
Уравнение (7.49) - одно из самых простых дифференциальных уравнений в частных производных. Несмотря на его элементарный вид, решение такого уравнения даже в простейшей ситуации есть весьма сложная задача.
Уравнение теплопроводности в трехмерном случае. Описанный выше вывод уравнения теплопроводности достаточно элементарен. Рассмотрим вывод уравнения теплопроводности в трехмерном случае, используя более общий аппарат математического анализа.
Рис. 7.31. Иллюстрация к выводу уравнения теплопроводности в трехмерном случае
Рассмотрим некоторое тело (V), ограниченное поверхностью (S) (рис. 7.31). Закон сохранения энергии должен выполняться для любой части тела (V). По этому закону скорость изменения энергии в теле равна потокуэнергии через его границу. Имеем для энергии в объеме V
где ? (
Поток энергии через границу тела S равен
Применяя к правой части теорему Остроградского- Гаусса, получаем
Поскольку это соотношение должно выполняться для любой части тела (V), то необходимо и достаточно, чтобы в любой точке
Учитывая, что плотность энергии ? (
где и = u (
Далее будет продолжено лишь рассмотрение задачи о теплопроводности в стержне.
Начальные и краевые условия.
Уравнения (7.49), (7.50) описывают процесс изменения температуры тела (перенос тепла) во времени и в пространстве. Ясно, что для отслеживания такого процесса надо знать распределение температуры в теле в некоторый начальный момент времени:
где f(x) - заданная функция. Кроме того, в тех местах, где возможен теплообмен с окружающей средой, надо знать условия этого теплообмена. Для стержня с теплоизолированной боковой поверхностью такими местами являются концы. Пусть длина стержня l; если один конец имеет координату x = 0, а.
другой - x
= l, то простейший вариант краевых условий - постоянная (но не обязательно одинаковая) температура на каждом конце стержня:
Нижеследующее утверждение физически очевидно, но его строгое математическое доказательство весьма непросто: дифференциальное уравнение (7.49) при начальном условии (7.51) и краевых условиях (7.52) имеет единственное решение.
Аналитические методы решения задачи одномерной теплопроводности существуют, но требуют значительной математической подготовки, к тому же решение обычно получается в виде ряда Фурье, и по его виду протекание процесса неочевидно. В двух- и трехмерном случаях аналитическое решение чаще всего получить не удается (по крайней мере, в практически полезном виде). Как и всюду в этой главе, ниже мы используем простейшие численные методы его решения. Вначале, однако, приведем графические результаты решений простейших задач (заимствованные из книги И.Г.Арамановича и В.И.Левина «Уравнения математической физики», Москва, 1969), способствующие пониманию рассматриваемой проблемы.
Пример 1. В конечном стержне (с теплоизолированной боковой поверхностью) оба торцевых сечения теплоизолированы, а начальная температура распределена по следующему закону:
Графики температуры построены в некоторые последовательные моменты времени, рис. 7.32. При любом t > 0 график симметричен относительно точки
Теплоизоляция концов стержня находит свое выражение в том, что кривые распределения температуры имеют горизонтальные касательные при
x = 0 и х = l. Из физических соображений ясно, что при t > ? u >uo/2.
Рис. 7.32. Графическая иллюстрация решения задачи из примера 1
Пример 2. В конечном стержне (с теплоизолированной боковой поверхностью) оба торцевых сечения теплоизолированы, а начальная температура распределена по следующему закону:
Здесь u0 - максимальное значение температуры.
В точках
Постоянная температура на торцах стержня - простейшее краевое условие. Возможна, однако, и ситуация, когда через торцы происходит теплообмен с окружающей средой. Этот теплообмен, как было установлено Ньютоном, удовлетворяет правилу: поток тепла через единицу поверхности в единицу времени пропорционален разности температур тела и окружающей среды: ?Q = h (u -
где и - температура конца стержня,
- коэффициент теплообмена. По определению h > 0, т.е. ?Q > 0 соответствует уходу тепла из стержня, ?Q < 0 - приходу из окружающей среды.
Рис. 7.33. Графическая иллюстрация решения задачи из примера 2
Поскольку поток тепла во внешнюю среду пропорционален градиенту изменения температуры на торце стержня, закон сохранения энергии принимает вид
(знак «минус» во второй формуле связан с соотношением направления потока и оси х), k - коэффициент теплопроводности.
Ниже приведен пример эволюции температуры в стержне, у которого один из концов теплоизолирован, а на другом - поддерживается постоянная температура.
Пример 3. В стержне (с теплоизолированной боковой поверхностью) левый конец теплоизолирован:
Рис. 7.34. Графическая иллюстрация решения задачи из примера 3
Методы конечных разностей в моделировании свойств сплошных сред. Покажем на примере уравнения теплопроводности наиболее распространенные методы численного интегрирования уравнении в частных производных. В их основе лежит прием дискретизации.
Покроем отрезок [а, b] одномерной сеткой (т.е. разобьем на n равных частей, рис. 7.35) с узлами в точках
Искомую функцию и(х)
будем аппроксимировать ее значениями в узлах сетки. Конечно, такое представление не дает полного описания, но в промежуточных точках, если сетка достаточно «мелкая», возможна интерполяция.
Рис. 7.35. Одномерная сетка
Остановимся на разностной аппроксимации производных. Производная дает информацию о локальном изменении функции в пространстве и, соответственно, связывает ее значения в соседних узлах сетки. Очевидная аппроксимация первой производной в точке х,
имеет вид
Для крайних точек, однако, такая аппроксимация невозможна, и простейший способ - ограничиться односторонними разностями:
Разумеется, (7.54) и (7.55) дают простейшие аппроксимации. Втягивая большое количество узлов, можно получить аппроксимации более высокого порядка, но часто бывает достаточно описанных выше. Аналогичная им аппроксимация вторых производных имеет вид
Что же касается методов интегрирования по времени, то это те же методы, что и для обыкновенных дифференциальных уравнений: Эйлера, Рунге - Кутта и т.д. Так как им тоже свойственна дискретизация, то возникает еще одна, временная сетка. При интегрировании уравнений по времени мы движемся по отдельным слоям, а в каждом слое определяем значение искомой функции на пространственной сетке. Если для интегрирования по времени используется метод Эйлера или другой одношаговый метод, то для работы со следующим временным слоем используются значения искомой функции из предыдущего слоя, для более сложных - из нескольких предыдущих слоев.
Далее будем индексы, соответствующие временной сетке, писать надстрочно (вверху), а пространственной - подстрочно (внизу).
Таким образом, для одномерного уравнения запись u
(k = 0, 1,...; i
= 1, 2, ..., n - 1) для внутренних узлов пространственной сетки; в силу начального условия (7.51)
где функция f(x) задана и определяет значение температуры при t = 0. Что касается значений u
Теперь остановимся на вопросе об устойчивости и эффективности обсуждаемого метода. Устойчивость понимается в том же смысле, что и для обыкновенных дифференциальных уравнений, но шансов получить неустойчивый метод здесь гораздо больше. Существуют разностные схемы абсолютно неустойчивые, абсолютно устойчивые и условно устойчивые. Первые при любых, сколь угодно малых, шагах так «раскачивают» начальную погрешность, что приводят к результатам, не имеющим ничего общего с реальностью. Вторые ни при каких шагах не «раскачиваются», хотя, конечно, чем меньше шаг, тем меньше разница между приближенным и точным решениями. Третьи устойчивы при одних комбинациях значений ?x и ?t и неустойчивы при других. Исследование, которого мы проводить не будем, показывает, что разностная схема (7.57) устойчива при
и неустойчива в противном случае.
Эффективность схемы можно представить лишь при сопоставления с другой схемой того же назначения. Прежде всего, под эффективностью понимают возможность относительно быстро получить решение с достаточной точностью. Иногда оказывается не менее важным объем оперативной памяти под массивы, хранение которых неизбежно в данном методе.
Схема (7.57) с точки зрения быстродействия малоэффективна, с точки зрения объема памяти - вполне удовлетворительна, так как, получив значения и
Получим более эффективный и устойчивый метод. Он аналогичен переходу от метода Эйлера к одному из вариантов метода Рунге - Кутта второго порядка (называемому иногда модифицированным методом Эйлера). Усредним пространственный член уравнения (7.49) по времени:
Это, безусловно, лучшая чем в (7.57) аппроксимация производной
Расплатой за эффективность является то, что (7.58) - неявная схема, т.е. не формула для непосредственного расчета, как (7.57), а система линейных алгебраических уравнений для величин u
Заметим, что (7.58) есть система специального вида - с трехдиагональной матрицей. В самом деле, если выписать первое, последнее и некоторое промежуточное ;'-е уравнения, перенося неизвестные в левые части, получим
Конечно, к таким системам можно применять стандартные методы решения систем линейных алгебраических уравнений, но для них существует и специализированный высокоэффективный метод, называемый «методом прогонки». За деталями отсылаем к учебникам по численным методам.
Пример.
Рассмотрим динамику изменения температуры в стержне длиной 4 м с теплоизолированными концами, температура на которых поддерживается постоянной и равна 3°С с начальным условием f(x) = -0,5x2 + 2x + 3. Коэффициент а в уравнении (7.49) примем равным 0,78 (выбор этот достаточно произволен).
Для демонстрации работы явной схемы (7.57) произведем расчеты по этой формуле на первом шаге.
Ограничимся пятью узлами на пространственной сетке. В начальный момент (t = 0) имеем u
4,5000, и
3,0000.
Из краевых условий получаем и
3,0000. Подставляя в формулу (7.57) соответствующие значения, получаем
аналогично получаем u
Таблица 7.7
Результаты моделирования процесса теплопроводности, полученные по неявной схеме (7.59)
 t |
0 |
1 |
2 |
3 |
4 |
|
0 |
3,000 |
4.500 |
5,000 |
4,500 |
3,000 |
|
1 |
3,000 |
4,000 |
4,428 |
4,000 |
3,000 |
|
2 |
3,000 |
3,688 |
3,975 |
3,688 |
3,000 |
|
3 |
3,000 |
3,476 |
3,669 |
3,476 |
3,000 |
|
4 |
3,000 |
3,325 |
3,461 |
3,325 |
3,000 |
|
5 |
3,000 |
3,225 |
3,316 |
3,225 |
3,000 |
|
6 |
3,000 |
3,154 |
3,218 |
3,154 |
3,000 |
|
7 |
3,000 |
3,106 |
3,150 |
3,106 |
3,000 |
|
8 |
3,000 |
3,073 |
3,103 |
3,073 |
3,000 |
|
9 |
3,000 |
3,050 |
3,071 |
3,050 |
3,000 |
|
10 |
3,000 |
3,034 |
3,049 |
3,034 |
3,000 |
Рис. 7.36. Графики зависимости температуры от координаты в разные моменты времени (сверху вниз t = 0, t = 2, t = 4, t = 6, t = 8), в начальный момент времени температура самая высокая, затем она постепенно выравнивается, и зависимости температуры от времени в разных точках стержня. Верхняя кривая соответствует x = 2; ниже - x = 1 и х = 3; прямая линия, совпадающая здесь с осью абсцисс, - значение температуры на концах стержня
Ясно, что по мере эволюции во времени температура стержня будет выравниваться и асимптотически стремиться к 3oС во всех точках.
