|
Лекция 11.
|
Рис. 11.1. Система взаимодействующих тел в задаче теплопроводности |
Как видно, в процессе жизни в системе изменяются (могут измениться) четыре показателя: температуры тел A, B, Tс, Tп. Значит, мы имеем дело с четырьмя переменными, зависящими от времени (поскольку переменные меняют свои значения со временем). Введем эти переменные: X1(t), X2(t), X3(t), X4(t).
Для построения математической модели данной системы отразим процесс теплопередачи в виде графа зависимостей (рис. 11.2).
Рис. 11.2. Граф зависимости переменных системы |
Если иметь в виду переменные X1(t), X2(t), X3(t), X4(t), то граф будет выглядеть так, как показано на рис. 11.3.
Рис. 11.3. Граф зависимости переменных модели |
Стрелка от A к B обозначает изменение температуры X2(t) объекта B под влиянием объекта A. Понятно, что ряд стрелок (например, от B к Tс, от A к Tп и др.) отсутствует, то есть нет влияния одних параметров на другие: тело B не в состоянии сколько-нибудь существенно нагреть открытую атмосферу, а тело A массивную и потенциально бесконечную опору. Строго говоря, такое влияние есть, но оно настолько ничтожно, что разумно им пренебречь.
Поскольку переменных четыре, то нам необходимо, как минимум, четыре закона, описывающих их изменение. В общем виде, учитывая, от каких переменных зависит каждый показатель, получим:
Стрелки, входящие в соответствующий кружок, указывают на количество влияющих параметров, а то, откуда они исходят, определяет конкретные названия переменных.
Для среды закон имеет вид: X3(t) = const, то есть, температура атмосферы Tс не зависит от остальных составляющих данной системы и, соответственно, не изменяется. Для опоры закон имеет вид: X4(t) = const, то есть, температура опоры Tп не зависит от остальных составляющих данной системы и, соответственно, не изменяется.
Система законов в первом приближении сформирована. Остается определить их конкретный вид: раскрыть, что из себя представляют значения выражений f1 и f2. Так как мы имеем дело с системой, зависящей от своего прошлого поведения на каждом последующем шаге, то мы применили для ее описания дифференциальные уравнения.
Основной динамический закон для описания изменения переменной (уравнение движения) имеет вид:
dX(t)/dt = w(x(t), y(t), z(t), ).
Физический смысл записи таков. Производная в левой части уравнения, по определению, показывает, насколько изменяется X с изменением времени t. В инженерии подобное изменение называется скоростью, темпом, тенденцией. Итак, чтобы записать закон изменения переменной в дифференциальных уравнениях, надо указать скорость изменения переменных.
Сначала рассмотрим первое уравнение:
dX1(t)/dt = f1(X1(t), X2(t), X3(t)).
Появление X1 в правой части означает, что скорость изменения температуры зависит от собственного состояния тела. Что такое f1? это функция, связывающая переменные X1, X2, X3 между собой. То есть, переменные соединены друг с другом знаками операций.
Взгляните на граф на рис. 11.3. Какие пары переменных взаимодействуют? Стрелки соединяют X1(t) с X2(t), X1(t) с X3(t), то есть имеет место два процесса, влияющих на скорость. Мы рассматриваем процессы теплообмена тел. Известно, что два процесса теплообмена независимы, то есть не управляют друг другом. Значит, результаты двух процессов можно складывать, они как бы накладываются друг на друга. Действительно, тепло, переданное от одного тела, складывается с теплом, переданным от другого. Таким образом, имеем:
dX1(t)/dt = g1(X1(t), X2(t)) + g2(X1(t), X3(t)).
Раскроем структуру оставшихся выражений g1 и g2. Очень удобно, что g1 никак не зависит от g2 и может рассматриваться отдельно. Такое разделение возможно, так как процессы g1 и g2 независимы. Процесс g1 идет независимо от того, идет или нет процесс g2. Независимость процессов и линейность (аддитивность) выражений понятия связанные. Итак: так как процессы g1 и g2 независимы, то забудем на некоторое время о g2.
Какой знак нужно поставить между X1(t) и X2(t) в выражении g1? Возможные варианты:
X1(t) + X2(t);
X2(t) X1(t);
X1(t) X2(t);
X1(t) · X2(t);
X1(t) / X2(t);
X2(t) / X1(t);
X1(t) ^ X2(t);
X2(t) ^ X1(t)
и далее более сложные, например, X12(t) · cos(X2(t))/exp(X1(t)). Исследователь начнет с наиболее простых выражений природа построена просто. И только если простейшие выражения не удовлетворяют исследователя, он переходит к более сложным вариантам описания.
То, что было только что сказано выше, возведено в системотехнике в ранг принципа: «не вводи сущностей без надобности» (принцип Оккама).
Итак, пусть: dX1(t)/dt = X1(t) + X2(t). Какие есть качественные варианты у этой физической системы?
Пусть теперь dX1(t)/dt = X2(t) X1(t). Какие есть качественные варианты у этой физической системы?
Других вариантов существования системы нет, и рассмотрение оканчивается.
Забыв на некоторое время о g1, так же можно рассмотреть и g2, что предлагается проделать читателю самостоятельно. В конечном итоге мы получим:
dX1(t)/dt = (X2(t) X1(t)) + (X3(t) X1(t)).
Далее. Так как, во-первых, у разных материалов разность температур влияет на скорость изменения температуры тела различным способом и, во-вторых, скорости двух процессов (у двух разных пар материалов) могут быть разными, то скорректируем модель, используя коэффициент теплопроводности, который играет роль усилителя (ослабителя) процессов. Это коэффициент влияния связи на объект. При K = 0 влияние отсутствует, связь отключается. При K = 0.0001 влияние слабое. При K = 1000 влияние связи огромно. Понятно, что коэффициент стоит при выражении процесса: K ? (X2(t) X1(t)), где «?» означает знак некоторой операции, а именно операции умножения. Эта операция дает зависимость одного члена от другого (в нашем случае K от X2(t) X1(t)).
При K = 0, какие бы значения ни принимало выражение X2(t) X1(t), результат K · (X2(t) X1(t)) дает 0. При K = 0.1 значение выражения принудительно ослабляется в 10 раз. То же самое верно и со стороны выражения X2(t) X1(t) это естественно, ведь общее выражение K · (X2(t) X1(t)) симметрично. Значит, мультипликативная связь моделирует свойство нелинейности и взаимозависимости процессов (один может сводить на нет действие другого). Сложение таким свойством не обладает.
В итоге модель имеет вид:
dX1(t)/dt = K21 · (X2(t) X1(t)) + K31 · (X3(t) X1(t)).
В конце следует проверить размерности уравнения; размерность левой части должна совпасть с размерностью правой. Напомним только, что производная имеет размерность показателя X, деленного на единицу времени.
Теперь мы в состоянии синтезировать аналогично второе уравнение (рекомендуем проверить правильность данного уравнения самостоятельно):
dX2(t)/dt = K12 · (X1(t) X2(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)).
Уравнение изменения температуры атмосферы: dX3(t)/dt = 0. То есть X3 = const (X3 не изменяется).
Уравнение изменения температуры опоры: dX4(t)/dt = 0. То есть X4 = const (X4 не изменяется).
Вся система уравнений в сборе имеет вид:
dX1(t)/dt = K21 · (X2(t) X1(t)) + K31 · (X3(t) X1(t));
dX2(t)/dt = K12 · (X1(t) X2(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t));
dX3(t)/dt = 0;
dX4(t)/dt = 0.
По физическим соображениям ясно: сколько тепла вытекает из A в B, столько же тепла поступает в В из А, то есть K21 = K12.
К записи общей модели остается добавить конкретные значения коэффициентов теплопроводности (K12 = K21, K31, K32, K42) и начальное состояние системы:
X1(0) = a,
X2(0) = b,
X3(0) = c,
X4(0) = d,
где a, b, c, d числа, указывающие температуру соответствующего объекта в момент времени t = 0.
|
Важное замечание. Проделанная выше работа по подтверждению гипотез, содержащихся в модели данного примера, естественно, не доказывает абсолютной правильности принятой модели. Проверки могут быть продолжены. Если при очередной проверке гипотеза будет отвергнута, то модель следует снова уточнить. Одной из дополнительных проверок может быть, например, проверка на открытость системы. Выполним такую проверку. Заметим, что выше мы получили открытую систему, то есть такую, чье суммарное тепло не постоянно, а может изменяться. Это видно из асимметрии стрелок на графе. Проверим этот факт математически, формально, для чего сложим левые части всех уравнений и, отдельно, правые части. Слева мы имеем следующее: dX1(t)/dt + dX2(t)/dt + dX3(t)/dt + dX4(t)/dt или d(X1(t) + X2(t) + X3(t) + X4(t))/dt или dXсистемы(t)/dt. В правой части мы имеем следующее: K21 · (X2(t) X1(t)) + K31 · (X3(t) X1(t)) + K12 · (X1(t) X2(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)) или (с учетом того, что K21 = K12) K31 · (X3(t) X1(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)). И, наконец, вместе: dXсистемы(t)/dt = K31 · (X3(t) X1(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)), из чего следует, что левая часть явно не равна нулю и есть утечка или приток тепла в систему извне.
В случае закрытой системы уравнения имели бы вид:
При сложении это дает: dXсистемы(t)/dt = 0 уравнение закрытой системы, в которую извне нечего не притекает и из которой ничего не истекает. Хотя внутри системы происходят изменения (перераспределение тепла), суммарная температура системы неизменна. Ни первая модель (модель открытой системы), ни вторая модель (модель закрытой системы) не являются «плохими». Просто достигнуты различные цели, даны различные представления исследователя о процессах, о системе в целом.
Динамическая система, которую мы рассматривали выше, это общий
случай. Весьма важно, что из нее всегда можно «даром» получить
статическую систему, для чего нужно потребовать:
Xi(t + Δt) = Xi(t),
в результате чего имеем:
dXi(t)/dt = 0.
Другими словами, запись
Xi(t + Δt) = Xi(t)
означает, что прошлое равно настоящему, то есть состояние системы не меняется.
Тогда от уравнений остается:
Это уравнения статики. После построения базовой системы исследователь может вводить дополнительные конструктивные элементы в систему тел. Добавим, например, нагреватель (см. рис. 11.4).
Граф изменится, так как на тело B будет действовать дополнительный объект нагреватель. Обозначим его температуру переменной X5(t) (см. рис. 11.5).
Изменения в записи будут касаться только второго уравнения, описывающего изменение переменной X2 во времени. Кроме того, нам понадобится уравнение нагревателя, так как в системе появилась дополнительная переменная X5(t) (вершина графа) и закон ее изменения должен быть определен. При записи уравнения следует различать, к какому типу относится источник тепла слабому или сильному. Рассмотрим сначала гипотезу о слабом источнике. Обратите внимание: если нагреватель будет греть сильнее тела B, то, конечно, тело B нагревается от нагревателя, но если нагреватель в какой-то момент окажется холоднее тела B, то тело B само отдает тепло нагревателю, как бы странно это, на первый взгляд, ни казалось. Так как такой источник энергии может не только нагревать другие тела, но и нагреваться от них сам, он называется слабым. То есть от услуг такого источника энергии можно отказаться, поскольку он находится в равном с другими телами положении, может влиять на них и может испытывать их влияние. На сильный источник другие тела воздействовать не могут, а сам он может на них воздействовать, НАВЯЗЫВАЯ ИМ СВОЮ ВОЛЮ БЕЗУСЛОВНО. Понятно, что в уравнении для тела B появится дополнительное слагаемое, аналогичное ранее рассмотренным: dX2(t)/dt = K12 · (X1(t) X2(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)) + K52 · (X5(t) X2(t)). Уравнение нагревателя может быть записано так: dX5(t)/dt = (K · U(t) X5(t))/L + K25 · (X2(t) X5(t)), где U(t) напряжение питания источника. Начальные условия X5(0) и функция работы источника U(t) должны быть заданы:
X5(0) = h,
Уравнение связывает причину (напряжение источника питания) и следствие (температура нагревателя). Свойства этого элемента таковы, что при выходе на рабочий режим нагреватель поддерживает постоянную температуру. При включении нагревателя рабочая температура устанавливается со временем, постепенно нарастая, при выключении источника нагреватель постепенно остывает. Очевидно, что процесс нагрева и остывания нагревателя инерционный. Фактически, для описания этих свойств нагревателя достаточно записать апериодический закон (который мы ранее уже обсуждали в лекции 04) изменения его температуры X5(t) по отношению ко входу U(t), что мы и сделали. В уравнении учтен коэффициент усиления K между входом U(t) и выходом X5(t), инерционность процесса нагрева L и коэффициент теплопередачи K25 между телами X2(t) и X5(t). Теперь рассмотрим гипотезу о сильном нагревателе. Это случай, когда нагреватель устроен таким образом, что излученная им энергия не может отразиться телом B, или если энергия от тела B не может попасть к нагревателю, или если тело B не может «отказаться» от услуг нагревателя (см. рис. 11.6).
На рис. 11.7 изображен граф зависимостей для данного случая.
Обратите внимание: неважно, нагрет ли нагреватель сильнее тела B или нет тело B будет получать от него энергию в любом случае, причем, в том объеме, в каком эту энергию будет выделять нагреватель. В уравнении тела B появится дополнительное слагаемое: dX2(t)/dt = K12 · (X1(t) X2(t)) + K32 · (X3(t) X2(t)) + K42 · (X4(t) X2(t)) + K52 · X5(t). Уравнение нагревателя может быть записано так: dX5(t)/dt = f(t). Начальные условия X5(0) и функция работы источника f(t) должны быть заданы: X5(0) = h. Но вернемся к расчету движения во времени динамической системы, для которой есть все необходимые данные, есть начальное состояние (Xi(0) = const). Можно по формулам (с использованием метода Эйлера, см. лекцию 10) вычислить скорость ее изменения и новое состояние: новое состояние := старое состояние + скорость · отрезок времени. Формально эта запись выглядит так: dx(t)/dt = f(x(t)) или (в дискретной форме) [x(t + Δt) x(t)]/Δt = f(x(t)) и, окончательно, x(t + Δt) = x(t) + f(x(t)) · Δt. Зная новое состояние и считая его уже достигнутым, то есть старым состоянием, используем ту же формулу (подставляя, конечно, в нее уже новые значения). И так продолжаем, прибавляя все новый отрезок времени, двигаясь во времени. Итак, мы перешли к вопросу о методе расчета дифференциальных моделей. Методом решения дифференциальных уравнений, в общем случае, является интегрирование их по независимой переменной времени t. Простейшим методом численного интегрирования дифференциальных уравнений является метод Эйлера (см. лекцию 10). Обратите внимание. При вычислении по методу Эйлера необходимо вычислять значения всех параметров системы параллельно, так как производная в каждый момент времени отдельного параметра зависит как от значения самого параметра в данный момент времени, так и от значения другого параметра в этот момент времени. Рассмотрим практически применение метода Эйлера для расчета процесса изменения температур тел системы. Зададим значения коэффициентов модели: K12 = K21 = 0.2, K31 = 0.1, K32 = 0.05, K42 = 0.1. Зададим начальные условия системы (в момент времени t = 0): X1(0) = 30° C, X2(0) = 70° C, X3(0) = 22° C, X4(0) = 15° C. Выбираем шаг моделирования Δt равный, например, 0.2 с. Примем конечное значение времени моделирования за Tk = 4 с.
Подставим значения коэффициентов:
и приступаем к моделированию. Процесс моделирования и численные значения отражены в табл. 11.1.
Во времени поведение системы будет выглядеть так, как показано на рис. 11.8.
Заметим, что на рис. 11.8 показана не вся траектория, а только ее часть. Рассчитайте всю траекторию, подумайте и ответьте на следующие вопросы:
В зависимости от реализации «машины», на которой автоматически будет имитироваться процесс, запись может выглядеть по-разному. Формальная математическая записьПриведем ниже вариант формальной математической записи для случая открытой системы тел со слабым нагревателем.
При реализации на математической машине среды Stratum эта запись будет дополнена элементами (функциями) ввода и вывода информации и командами управления (например, «стоп»), а часть элементов будет спрятана (например, при задании начальных условий). Иногда (когда это удобно) начальные условия могут быть помещены внутрь уравнений. Для этого в следующей записи мы воспользуемся дельта-функцией Дирака (см. лекцию 31).
Напомним, что osc2d(T, X1) функция, рисующая на экране точку с координатами T и X1 для каждого момента времени. Выше данная функция приведена условно: в действительности надо использовать несколько функций для настройки окон, масштаба изображения, цвета, толщины, стиля линии и тому подобного. Подробно реализацию двухмерного осциллографа в среде «Stratum-2000» вы можете посмотреть в тексте имиджа OSCSpace2D. Приведенный выше вариант формальной математической записи подразумевает, что весь код заключен в рамках одного элемента (имиджа, если использовать термины среды «Stratum-2000»). Будет гораздо более наглядно, если распределить код по отдельным элементам, связав их между собой связями (см. рис. 11.9) в этом случае структура проекта останется обозримой; кроме того, если потребуется изменить что-то в одном из блоков, то такой подход позволит не менять код в остальных блоках. Еще один плюс такого подхода отдельные «кирпичики» схемы могут быть применены неоднократно в готовом виде в виде копии имиджа.
Алгоритмическая реализацияВ случае реализации на алгоритмической машине следует воспользоваться алфавитом алгоритмов: блоками начала и конца алгоритма, ввода-вывода информации, задания начальных условий, вычислительными блоками, блоками условия, конструкциями цикла. На рис. 11.10 представлена блок-схема алгоритма, имитирующего теплообмен между телами.
Заметьте, что для реализации нам понадобились следующие инструменты:
В теле цикла присутствуют:
В той или иной мере вышеназванные блоки обязательно присутствуют в каждом алгоритме данного типа. |
Лекция 10. Численные методы интегрирования | Лекция 12. Уравнения высших порядков | ||||||||||||||||
|