Напоминаем, где находится VBA.
Литература: [0 – 3].
Основные понятия и типы моделей, способы их наглядного представления
Процессов теплопередачи.
Цель:набирать объем понятий для моделирования теплопередачи (в частности) и сплошных сред (вообще).
 |
Задача 1. Термостатика: моделирование установившегося распределении температур. Определите потери тепла через наружную стену здания при различных комбинациях стройматериалов (рис. 1), теплопотери трубы с горячей водой (рис. 2) и теплопотери через закругленный угол здания (рис. 3).
Для трехслойной стены (рис. 1, в) подберите толщину теплоизоляционного слоя так, чтобы теплопотери через нее были на 20 % ниже, чем через кирпичную стену (рис. 1, а).
Подсказкик рис. 1, в.Составьте систему уравнений, состоящую из законов теплопроводности для слоев (рис. 2, 3):
px0 = 0.07(T1 – T0)/0.25, px1 = 0.7(T2 – T1)/0.5, px2 = 0.07(T3 – T2)/0.25 (рис. 2, закон теплопередачи Фурье[2]);
px0 = px1 , px1 = px2 (рис. 3, закон сохранения энергии) ;
T0 = 20°, T3 = –20° (рис. 1 – 3, граничные условия).
Подсчитайте количество неизвестных и уравнений, должно совпасть.
Используйте подручное вычислительное средство – электронные таблицы MS Excel. Разместите на листе Excel имена неизвестных, номера уравнений, коэффициенты при неизвестных и правые части. Используйте функции вычисления обратной матрицы МОБР и умножения матриц МУМН для решения системы уравнений.
Эти уравнения можно решить последовательной подстановкой. Но это – не наше барское дело. Для этого у нас электронный аппарат – компьютер, надо только передать ему всю информацию. Используйте матричные функции Excel: мобр (вычисление обратной матрицы) и мумн (умножение матриц).
Вспомните: информация о системе линейных алгебраических уравнений состоит из двух таблиц чисел: матрицы коэффициентов и столбца (или нескольких столбцов) правых частей.
При размещении табличной информации на листе Excel целесообразно в верхней строке и левом столбце подписать заголовки: имена неизвестных и номера уравнений (см. табл. 1). Это нужно, чтобы без ошибок заполнить таблицы числами или формулами.
Для размещения обратной матрицы выделяют прямоугольную область из пустых ячеек n´n (см. табл. 2), и нее вставляют функцию вычисления обратной матрицы. Эта функция показывает диалоговое окно «Аргументы». В нем нужно указать адреса левой верхней и правой нижней ячеек (например, B3 : H9). Или выделить (выбрать) этот диапазон ячеек с помощью мышки (если он помещается на экране).
Следующий важный момент: кнопку OK в диалоговом окне не нажимают[3], а держат левой рукой клавиши Ctrl и Shift, а правой рукой нажимают и отпускают клавишу Enter (только после этого отпускают Ctrl и Shift)[4].
| Таблица 1 Система линейных алгебраических уравнений: матрица коэффициентов при неизвестных и столбец правых частей | |||||||||
| T0 | T1 | T2 | T3 | p0 | p1 | p2 | b | ||
| у1 | -0,28 | 0,28 | |||||||
| у2 | -1,4 | 1,4 | |||||||
| у3 | -0,28 | 0,28 | |||||||
| у4 | -1 | ||||||||
| у5 | -1 | ||||||||
| у6 | |||||||||
| у7 | -20 |
Получив обратную матрицу, таким же путем в свободный столбец из n ячеек (или в несколько столбцов, по количеству столбцов правых частей) вставляют функцию умножения матриц, указывают в появившемся диалоговом окне ее аргументы. Затем нажимают Ctrl+Shift – Enter и получают столбец неизвестных в том порядке, в котором их имена стоят в строке заголовков матрицы коэффициентов.
| Таблица 2 Обратная матрица и столбец неизвестных | x | x | ||||||
| T0 | ||||||||
| 1,948052 | -0,32468 | -1,62338 | -1,94805 | -1,62338 | 0,545455 | 0,454545 | 1,818182 | T1 |
| 1,623377 | 0,324675 | -1,94805 | -1,62338 | -1,94805 | 0,454545 | 0,545455 | -1,81818 | T2 |
| -20 | T3 | |||||||
| 0,454545 | 0,090909 | 0,454545 | 0,545455 | 0,454545 | 0,127273 | -0,12727 | 5,090909 | p0 |
| 0,454545 | 0,090909 | 0,454545 | -0,45455 | 0,454545 | 0,127273 | -0,12727 | 5,090909 | p1 |
| 0,454545 | 0,090909 | 0,454545 | -0,45455 | -0,54545 | 0,127273 | -0,12727 | 5,090909 | p2 |
 |
Изобразите наглядно результаты в виде эпюр (графиков) и изополей температуры и теплового потока (рис. 4).
Задача 2 . Термостатика: моделирование установившегося распределении температур при выделении тепла внутри объема.
Определите, до какой максимальной температуры нагреется свежезалитый бетонный блок за счет тепловыделения при гидратации цемента (рис. 5).
Подсказки. Очевидно, со временем в центре блока должна установиться максимальная температура Tmax и нулевая скорость теплопередачи p (тепловой поток), а на торцах – максимальная скорость теплопередачи pmax, так как при установившемся («стационарном») состоянии вся выделяемая тепловая энергия передается во внешнюю среду через торцы блока, вправо и влево одинаково, симметрично. Значит, распределение T должно быть нелинейным (рис. 5, г), и для определения его и тепловых потоков нужно использовать основной метод математического анализа – дифференцирование, т.е. разделение на части, моделирование по частям.
Используйте симметрию для уменьшения количества неизвестных и уравнений (рис. 5, а, в). Простейшая техника учета симметрии такова. Неизвестные (искомые) значения температур в симметричных точках обозначают одинаковыми именами (на рис. 5, в: T0 и T0 , T1 и T1 ). Все уравнения (закон теплопередачи Фурье и баланс тепловых потоков) составляют только для одной симметричной половины блока (например, для левой). В некоторые из этих уравнений попадают значения T из правой
 |
половины, но это не новые неизвестные, количество неизвестных и уравнений совпадает.
Но обычно предпочитают усовершенствованный вариант учета симметрии, технически более простой и универсальный. Так как температуры слева и справа от плоскости симметрии одинаковы, то тепловая энергия через эту плоскость не идет. Значит, условие симметрии равносильно условию теплоизоляции pS=0 (рис. 6, а).
А. Для первого самого грубого приближения представьте всю левую половину блока как одну крупную часть (рис. 7, а). Для нее составьте зависимость теплового потока p0.5от разности соседних температур T1 и T0 (закон теплопередачи Фурье, рис. 6, а), как в задании 1. Рис. 7, б иллюстрирует, каким образом получается, что для плавно меняющихся функций разность двух значений, деленная на расстояние между аргументами, дает с повышенной точностью значение производной (т.е. главного коэффициента линейного приближения) в центре интервала. Нужно либо рисовать чертежи и писать формулы на бумаге, либо делать это на компьютере [5].
К уравнениям теплопередачи надо добавить условия баланса тепловых потоков – закон сохранения энергии (рис. 7, б или в). Вспомните, что мы строим модель установившегося распределения температур, т.е. того состояния тепловых потоков, когда температура в каждом малом объеме уже перестала изменяться, наступило равновесие (баланс) между выделяющейся в этом объеме тепловой энергией и выходящей из него через его поверхность.
В нашем приближенном решении наиболее подходящий объем для составления такого баланса ограничивается сечением 0.5, в котором мы имеем наиболее точное значение p = p0.5 (пока неизвестное, выраженное через неизвестную температуру T1) и сечением 1, в котором p = pS известен, равен нулю.
Для составления баланса обычно приравнивают нулю сумму всех энергий, добавляющихся в выбранный объем за счет тепловыделения (W×DV1) и через поверхности p0.5×A – pS×A.[6] В результате уравнение установившегося баланса потоков энергий имеет вид
p0.5×A – pS×A + W×DV1 = 0. (1)
Полученную систему двух уравнений с двумя неизвестными можно решить в электронных таблицах Excel, используя матричные функции мобр и мумн, или вручную.
 |
Б.Для второго, более точного приближения разделите левую симметричную половину блока на две части (рис. 8, а), для каждой из них составьте условие теплопередачи Фурье. Для трех блоков, ограниченных поверхностями с неизвестными потоками энергии p0, p0.5, p1.5 и известным потоком p2=0, составьте условия баланса этих потоков и энергии тепловыделения в бетоне. Решите полученную систему уравнений с помощью матричных функций MS Excel и оформите отчет о лабораторной работе средствами MS Word.
Оформление.Представьте в документе Word исходные текст задания, исходные данные, процесс моделирования и результаты. Результаты изобразите наглядно в виде эпюр температурного поля и тепловых потоков. [7]
Результат работы.Представления о температуре T(x) как мере неравномерности нагрева частей тела, неравномерности распределения тепловой энергии. Эту неравномерность в окрестности каждой точки тела измеряют производными от температуры по разным направлениям. Производные в данной точке по всем направлениям определяются вектором производной g = dT/dx, т.е. на плоскости парой чисел (gx = dT/dx, gy = dT/dy), а в 3D пространстве тройкой чисел – проекций на оси декартовых координат. Этот вектор называют также градиентом функции T(x, y, z). От неравномерности распределения температуры T, т.е. от ее градиента, зависит движение энергии, измеряемое векторами интенсивностей (или плотностей) потоков энергии p(x, y, z). Плотность потока – это количество (энергии или вещества), проходящее за единицу времени через единицу площади поверхности, перпендикулярной исследуемому направлению потока. Так что px – это поток энергии в направлении x через поверность, перпендикулярную оси x, а py – через поверхность, перпендикулярную оси y.
Задача 3 . Термодинамика –моделирование неустановившегося поля температур, расчет его изменения. Оформление, наглядное представление процесса.
1. Моделирование и расчет. Успеет ли бетонный блок (рис. 5 к задаче 2) разогреться до своих примерно 200 или 300 градусов за те полторы – две недели, когда идет процесс гидратации цемента?В каком месте и какой величины будет максимальная температура через две недели после изготовления блока?
За какое время температура блока приблизится к установившейся температуре, найденной в лабораторной работе 1?
Другие примеры термодинамики: глубина промерзания грунта;
толщина льда в реке: дорога в Норильск, Дудинку, Туру;
огнезащита металлических конструкций
Подсказки. Предыдущие задачи относились к термостатике – расчету установившегося распределения температур (при наличии источников и стоков тепловой энергии). Это – относительно редкий частный случай, в обычных условиях температура тел постоянно меняется: увеличивается при добавлении тепловой энергии и уменьшается при убывании.
Используйте рис. 8, рис. 9 и уравнение (1) для составления уравнений, определяющих скорость изменения количества тепловой энергии в каждом малом объеме блока и соответствующего изменения температуры этого объема. Вспомните, что левая часть в (1) представляет собой сумму всех тепловых энергий, добавляемых в секунду в объем V1 (т.е. суммарная мощность добавления тепла, Вт).
 |
Например, в начальный момент времени t0 для левого участка блока (вокруг температуры T0, рис. 7, б, в) величина
(Вт) (здесь А – площадь сечения блока, А = 1 м2, рис. 5) (2)
означает добавление (приращение) тепловой энергии E в этот участок за секунду, т.е. скорость ее изменения (мощность добавления энергии). Чтобы смоделировать установившееся, постоянное тепловое состояние блока, мы эту величину приравнивали к нулю. Она (по условию задачи) равна нулю и в процессе изменения температурного поля – по условию задачи (дует ветерок, уносит тепло и поддерживает постоянное значение T0).
Для моделирования процесса изменения поля температур (а с ним и тепловых потоков) во втором узле моделируемого блока будем использовать скорость изменения E (например, в объеме вокруг точки T1) для получения значений E (а с ними и значений T1) в следующий момент времени:
, (3)
. (4)
Так как температура – мера тепловой энергии в единице объема вещества (с учетом его физических свойств – коэффициента теплоемкости с и плотности массы r), то по величине 
можно определить добавку DT1 температуры в этом объеме за тот же промежуток времени Dt:
(где m – масса, m = rV ) (5)
(по смыслу коэффициента теплоемкости материала c, Дж/кг), или
или 
, где DE1 вычисляется по (3). (6)
Дополнительные разъяснения. Уравнения термодинамики составляют из тех же понятий, что и уравнения термостатики. Кратко напомним их.
Температура T – мера уровня энергии, влияющего на теплообмен с соседними телами, на передачу ее от данного тела к смежному или обратно. Теплообмен зависит не от величины температуры, а от разностей температур в различных точках сплошной среды (сравните с понятиями электрического потенциала и разностей потенциалов). Разность температур (как и электрических потенциалов) между точками сплошной среды величина направленная (векторная). В моделировании не используют любые пары точек (это слишком много и сложно), а применяют дифференцирование с линейным приближением функции распределения температур T(x,y) в окрестностях точек. Вместо разности температур используют удельную разность (на единицу расстояния), т.е. производную 
. Здесь ¶x – величина направленная (векторная), поэтому и производная g– вектор, 
, ее называют градиентом[8]функции T(x,y).
Теплообмен, или передача тепловой энергии от боле нагретых частей (с более высокой температурой) к менее нагретым частям (т.е. в сторону убывания температуры), происходит через поверхности, разделяющие эти части. Теплообмен характеризуется понятием интенсивности (плотности) потока тепловой энергии: сколько энергии переходит из одной части в другую через единицу площади в единицу времени, Дж/(м2×сек) или Вт/м2. Интенсивность теплообмена p зависит от градиента температуры g и от свойства материала – коэффициента теплопроводности a: p=–a×g. Знак «минус» показывает, что тепло переходит в направлении убывания температуры (а градиент функции направлен в сторону ее наискорейшего возрастания).
Если материал анизотропный, то его теплопроводность представляют набором (матрицей) коэффициентов a. В этом случае направления векторов p и g чаще всего разные (немного различаются).
Можно было бы разделить обе части (5) на AdxDt и записать (5) в привычной форме дифференциального уравнения с частными производными (по x и по t), но для нашей первоначальной цели в этом нет никакой необходимости. Это даже помешает учесть особенность крайних участков блока – меньшую длину dx/2.
Уравнение (5) составлено для участка длиной dx, выделенного из всего моделируемого блока и расположенного вокруг ключевой точки с неизвестной температурой T1. Такие же уравнения, как (5), составляют для каждого блока длиной dx вокруг остальных неизвестных температур и из них вычисляют температуры T в каждом участке в следующий момент. (Длина последнего участка dx/2). Такие же вычисления проделывают для следующих моментов, и т.д.
Таким образом, в задачах динамики не решают системы уравнений[9]. Использование математической модели сводится к вычислению (из уравнений вида (3)) состояния объекта моделирования в следующий момент времени на основе уже имеющегося состояния в данный момент.
Для реализации модели этого процесса можно использовать формулы в ячейках таблицы Excel. Для этого разместите в электронной таблице Eхcel исходные данные. В одной строке поместите начальные значения времени t и температуры T в выбранных точках блока (узлах сетки). В следующей строке – формулу t+Dt и формулы (3) для вычисления следующего значения T через предыдущие. Скопируйте эту строку вниз на нужное количество суток и анализируйте результаты.
Однако это малоэффективный способ, он хорош только для предварительных проверок. Например, для достижения нужной точности модели может понадобиться выполнить миллион шагов по времени. В электронных таблицах не предусмотрено столько места, максимальный размер таблицы 64 тыс. строк и столбцов. Кроме того, как вы будете разбираться в таком огромном количестве данных? Только на прокручивание таблицы будут уходить сутки и недели. А для практического анализа процесса и применения в проекте нужно всего от 100 до 1000 строк.
Эффективный и универсальный путь – составление программы, которая будет вычислять значения T и p в последовательные моменты времени с любым, как угодно малым, шагом, а запоминать и отображать в таблице для осмысления, анализа и проверки только ограниченное, всегда одинаковое количество моментов t (например, 100).
Важное преимущество программы: все действия и формулы перед глазами, легче продумывать и уточнять порядок действий, сравнивать и проверять формулы.
Такую программу просто, легко, и удобно составлять в прикладной среде программирования VBA (Visual Basic for Application), связанной с электронными таблицами Excel, используя эти таблицы для размещения исходных данных и результатов моделирования, а также для обработки результатов (наглядного изображения в виде графиков и диаграмм, для простых вычислений по мере необходимости).
Для получения достоверного результата здесь требуется настройка величин интервалов dx, dt и их согласование: dt должен быть достаточно малым, чтобы линеаризованная по (3) траектория изменения T не слишком отклонялась от реальной. При уменьшении dx уменьшают и dt.
Из опыта авторов. Моделируем процесс разогрева сечения стальной конструкции с огнезащитой (рис. 7б, 1). Используем очень популярную программу моделирования сплошных сред. В ней по умолчанию автоматический выбор шага Dt. Результат в десять раз отличается от ожидаемого: не справилась программа со взятым на себя обязательством. Ничего особенного, в такой огромной работе не могут отсутствовать мелкие неувязки.
Реализуем упрощенную грубо приближенную модель (рис. 7б, 2), подбираем Dt. Результат близкий к ожидаемому. Настраиваем популярную программу на поиск шага Dt и получаем достоверный результат.
Назначьте максимальный разумный (по Вашему мнению) интервал dx. Начальный интервал dt выберите исходя из ожидаемого срока разогрева (например, t =1 месяц). Можно попробовать Dt = 1 сутки или Dt = 1 час. Выполните расчет и запомните результаты (например, скопируйте нужные колонки в другое место таблицы). Теперь задайте вдвое меньший интервал Dt, выполните расчет и сравните эти два результата. Если в конце процесса разница не превышает 5%, то это достаточная точность для наших целей, если нет – запомните последний результат и снова уменьшите Dt. Повторяйте, пока не получите достаточно малую разницу.
Чтобы легче сравнивать два процесса с разными Dt, желательно выводить в таблицу не все результаты расчета, а в одинаковые моменты для всех сравниваемых процессов. Для этого в исходных данных назначайте не только шаг Dt для расчета, но и шаг Dt_out для вывода. Вариант реализации этого усовершенствования представлен в табл. 2.
Используйте образец ниже: максимально упрощенный (табл. 1) или более длинный, но удобный для анализа и отображения результатов (табл. 2).
| Таблица 1 (см. также файл Образец_1_ПроцессаБэйсик.xls) Упрощенный образец программы для моделирования процесса изменения температур и тепловых потоков со временем на VBA (заготовка) |
| Private Sub CommandButton1_Click() Dim time As Double ' Const ro = 2400: c = 840: dtime = 60# * 60 * 24: alfa = 1.5: dx = 1.5 Const W = 100 Dim V, T, dT, p, Headers As Variant ', T As variant, dT As variant V = Array(0.75, 1.5, 0.75): T = Array(-20#, 0#, 0#): dT = Array(0#, 0#, 0#) p = Array(0#, 0#, 0#) Headers = Array("time,сутки", "", "T(0),грC", "T(1)", "T(2)", "", "p(0),Вт/м2", "p(1)", "p(2)") For i = 0 To 8: Лист1.Cells(2, i + 1) = Headers(i): Next i ' Вывод заголовков в ячейки таблицы Лист1 (в строку 2) Лист1.Cells(3, 1) = time / (24# * 3600#) ' Вывод начальных значений time, T, p в строку 3 For i = 0 To 2: Лист1.Cells(3, 3 + i) = T(i): Лист1.Cells(3, 7 + i) = p(i): Next i For i = 1 To 360 ' Цикл по промежуткам времени p(1) = -alfa * (T(1) - T(0)) / dx ' Предварительные вычисления : p(1) из закона теплопередачи Фурье p(2) = -alfa * (T(2) - T(1)) / dx ' p(2) из закона теплопередачи Фурье p(0) = p(1) - W * V(0) ' p(0) из баланса энергии в объеме V(0): когда уже определен p(1) c1 = (p(1) * 1 - p(2) * 1 + W * V(1)) * dtime c2 = V(1) * ro * c ' Коэффициенты основного уравнения - баланса энергии (для V(1)) dT(1) = c1 / c2 ' Определение dT(1) из этого основного уравнения c1 = (p(2) * 1 - 0 + W * V(2)) * dtime c2 = V(2) * ro * c ' Коэффициенты основного уравнения - баланса энергии (для V(2)) dT(2) = c1 / c2 ' Определение dT(2) из этого основного уравнения For j = 0 To 2: T(j) = T(j) + dT(j): Next j ' Добавление приращений температур в узлах сетки КЭ ' за время dtime к предыдущим температурам time = time + dtime Лист1.Cells(3 + i, 1) = time / (24# * 3600#) ' Вывод начальных значений time, T, p в строку 3+i For j = 0 To 2: Лист1.Cells(3 + i, 3 + j) = T(j): Лист1.Cells(3 + i, 7 + j) = p(j): Next j Next i End Sub |
Образцы составлены для минимальной сетки неизвестных температур (рис. 6). Скопируйте текст выбранного вами образца и проверьте работоспособность в среде Excel VBA[10]. Порядок действий может быть таким (для Excel 2007).
1) в MS Excel 2007. В главном меню активизируйте вкладку «Разработчик» и на ней нажмите кнопку “Visual Basic” для активизации окна программирования на Бэйсике.
Если в главном меню Excel нет вкладки «Разработчик», ее нужно активизировать, нажав кнопку “Office”, в появившемся меню (внизу) кнопку “Параметры Excel” и в следующем меню поставив «галочку» у надписи «Показывать вкладку “Разработцик” на ленте».
2) в MS Excel 97 – 2003.Цепочка пунктов главного меню: Сервис – Макрос – Редактор Visual Basic.
В главном меню окна «Редактор Visual Basic» выберите пункт “Insert” (Вставить) и в нем подпункт “User Form”. Появится заготовка окна будущей программы и рядом с ней “Toolbox” («Ящик с инструментами») Выберите в Toolbox (щелчком мышки) кнопку CommandButton и щелкните мышкой на форме, чтобы поставить на нее выбранную кнопку. Это – управляющий визуальный компонент операционной системы Windows для запуска процедуры решения задачи термостатики. Используя окно свойств программных объектов, переименуйте свойство кнопки Name (имя) в . Теперь сделайте двойной щелчок на этой кнопке, чтобы привязать событие контроллера прерываний TermoStatics к запуску подпрограммы TermoStatics_Click().
Поставьте на форму другую кнопку, дайте ей имя TermoDynamics и двойным щелчком мыши привяжите событие к запуску подпрограммы TermoDynamics_Click().
Запустите программу (нажав кнопку Run в главном меню VBA). Проверьте работоспособность подпрограмм термостатики и термодинамики, запустив их по очереди с помощью соответствующих кнопок и оценив результаты, помещенные ими на Лист1 и Лист2 табиц Excel.
Теперь сделайте еще одну копию и переделайте ее в программу расчета процесса для уточненной модели – с более густой сеткой неизвестных температур и тепловых потоков. Для этого нужно увеличить количество неизвестных и уравнений в глобальных определениях структур данных (в верхней части текста). В подпрограмме термостатики нужно заменить коэффициенты уравнений, увеличив их количество. В подпрограмме термодинамики нужно изменить операторы присваивания, вычисляющие необходимые элементы уравнений в нужном порядке.
Этот же образец понадобится в лабораторных работах 3 и 4, а также в контрольных заданиях 2 и 3.