Практика 6. Составление уравнений динамической системы (линейной и нелинейной) и расчет ее поведения методом Эйлера

Задание 1

Выполним построение модели динамической системы в виде дифференциальных уравнений и расчет ее методом Эйлера на примере.

Пример. Пусть исследуется система двух материальных тел A и B с различными теплофизическими свойствами. Система контактирует с опорой с температурой T_п и помещена во внешнюю среду с температурой T_c. Интересует протекание процесса изменения температур тел.

Рисунок - Система взаимодействующих тел в задаче теплопроводности

Как видно, в процессе жизни в системе изменяются (могут измениться) четыре показателя: температуры тел A, B, T_с, T_п. Значит, мы имеем дело с четырьмя переменными, зависящими от времени (поскольку переменные меняют свои значения со временем). Введем эти переменные: X₁(t), X₂(t), X₃(t), X₄(t).

Для построения математической модели данной системы отразим процесс теплопередачи в виде графа зависимостей.

Рисунок - Граф зависимости переменных системы

Если иметь в виду переменные X₁(t), X₂(t), X₃(t), X₄(t), то граф будет выглядеть так, как показано на рисунке.

Рисунок - Граф зависимости переменных модели

Стрелка от A к B обозначает изменение температуры X₂(t) объекта B под влиянием объекта A. Понятно, что ряд стрелок (например, от B к T_с, от A к T_п и др.) отсутствует, то есть нет влияния одних параметров на другие: тело B не в состоянии сколько-нибудь существенно нагреть открытую атмосферу, а тело A — массивную и потенциально бесконечную опору. Строго говоря, такое влияние есть, но оно настолько ничтожно, что разумно им пренебречь.

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

• для тела A имеем зависимость температуры X₁(t) от температуры тела B и температуры атмосферы T_с: dX₁(t)/dt = f₁(X₂(t), X₃(t));
• для тела B имеем зависимость температуры X₂(t) от температуры тела A, температуры атмосферы T_с и опоры T_п: dX₂(t)/dt = f₂(X₁(t), X₃(t), X₄(t)).

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

Для среды закон имеет вид: X₃(t) = const, то есть, температура атмосферы T_с не зависит от остальных составляющих данной системы и, соответственно, не изменяется. Для опоры закон имеет вид: X₄(t) = const, то есть, температура опоры T_п не зависит от остальных составляющих данной системы и, соответственно, не изменяется.

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

Основной динамический закон для описания изменения переменной (уравнение движения) имеет вид:

dX(t)/dt = w(x(t), y(t), z(t), …).

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

Сначала рассмотрим первое уравнение:

dX₁(t)/dt = f₁(X₁(t), X₂(t), X₃(t)).

Появление X₁ в правой части означает, что скорость изменения температуры зависит от собственного состояния тела. Что такое f₁? — это функция, связывающая переменные X₁, X₂, X₃ между собой. То есть, переменные соединены друг с другом знаками операций.

Взгляните на граф на рисунке. Какие пары переменных взаимодействуют? Стрелки соединяют X₁(t) с X₂(t), X₁(t) с X₃(t), то есть имеет место два процесса, влияющих на скорость. Мы рассматриваем процессы теплообмена тел. Известно, что два процесса теплообмена независимы, то есть не управляют друг другом. Значит, результаты двух процессов можно складывать, они как бы накладываются друг на друга. Действительно, тепло, переданное от одного тела, складывается с теплом, переданным от другого. Таким образом, имеем:

dX₁(t)/dt = g₁(X₁(t), X₂(t)) + g₂(X₁(t), X₃(t)).

Какой знак нужно поставить между X₁(t) и X₂(t) в выражении g₁? Возможные варианты:

X₁(t) + X₂(t);

X₂(t) - X₁(t);

X₁(t) - X₂(t);

X₁(t) • X₂(t);

X₁(t) / X₂(t);

X₂(t) / X₁(t);

X₁(t) ^ X₂(t);

X₂(t) ^ X₁(t);

и далее более сложные, например, X₁²(t) • cos(X₂(t))/exp(X₁(t)). Исследователь начнет с наиболее простых выражений — природа построена просто. И только если простейшие выражения не удовлетворяют исследователя, он переходит к более сложным вариантам описания.

То, что было только что сказано выше, возведено в системотехнике в ранг принципа: «не вводи сущностей без надобности» (принцип Оккама).

Итак, пусть: dX₁(t)/dt = X₁(t) + X₂(t). Какие есть качественные варианты у этой физической системы?

• X₁(t) > X₂(t). Тело A теплее тела B. Теплопоток при контакте двух тел направлен от A к B. Тело A отдает тепло телу B. То есть в процессе контакта значение X₁(t) падает — уменьшается. (Нас интересует будущее именно X₁(t), а не X₂(t) — см. уравнение: dX₁(t)/dt). Посмотрим, так ли это в уравнении-гипотезе dX₁(t)/dt = X₁(t) + X₂(t)? Сумма X₁(t) + X₂(t) может принимать как положительные, так и отрицательные значения, следовательно, значение dX₁(t)/dt также может быть как положительным, так и отрицательным, а это, в свою очередь, значит, что X₁(t) то растет, то падает. Но это противоречит физической картине, рассмотренной чуть выше: мы заключили, что при условии X₁(t) > X₂(t) X₁(t) может только лишь уменьшаться. Поэтому вариант гипотезы dX₁(t)/dt = X₁(t) + X₂(t) неприемлем и надо пробовать другой.

Пусть теперь dX₁(t)/dt = X₂(t) – X₁(t). Какие есть качественные варианты у этой физической системы?

• X₁(t) > X₂(t). Тело A теплее тела B. Теплопоток при контакте двух тел направлен от A к B. Тело A отдает тепло телу B. То есть в процессе контакта значение X₁(t) падает — уменьшается. Посмотрим, так ли это в уравнении? X₁(t) > X₂(t), то есть X₂(t) – X₁(t) < 0, значит, dX₁(t)/dt < 0, следовательно, X₁(t) падает. Вывод не противоречит физической картине. Значит, пока данный вариант приемлем и надо проверить его на остальных качественных ситуациях.
• X₁(t) < X₂(t). Тело A холоднее тела B. Теплопоток при контакте двух тел направлен от B к A. То есть в процессе контакта значение X₁(t) растет — увеличивается. Посмотрим, так ли это в уравнении? X₁(t) < X₂(t), то есть X₂(t) – X₁(t) > 0, значит, dX₁(t)/dt > 0, следовательно, X₁(t) растет. Вывод не противоречит физической картине. Значит, пока данный вариант приемлем и надо проверять его далее.
• X₁(t) = X₂(t). Температура тела A равна температуре тела B. Теплопоток при контакте двух тел равен нулю. То есть значение X₁(t) не изменяется — тело A не отдает и не принимает тепло. Посмотрим, так ли это в уравнении? X₁(t) = X₂(t), значит, X₂(t) – X₁(t) = 0, значит, dX₁(t)/dt = 0, значит, X₁(t) не изменяется. Вывод не противоречит физической картине. Значит, данный вариант принимается, так как он правильно (пока только качественно!) отражает физическую картину во всех случаях.

Других вариантов существования системы нет, и рассмотрение оканчивается.

Забыв на некоторое время о g₁, так же можно рассмотреть и g₂, что предлагается проделать читателю самостоятельно. В конечном итоге мы получим:

dX₁(t)/dt = (X₂(t) – X₁(t)) + (X₃(t) – X₁(t)).

Далее. Так как, во-первых, у разных материалов разность температур влияет на скорость изменения температуры тела различным способом и, во-вторых, скорости двух процессов (у двух разных пар материалов) могут быть разными, то скорректируем модель, используя коэффициент теплопроводности, который играет роль усилителя (ослабителя) процессов. Это коэффициент влияния связи на объект. При K = 0 влияние отсутствует, связь отключается. При K = 0.0001 влияние слабое. При K = 1000 влияние связи огромно. Понятно, что коэффициент стоит при выражении процесса: K ? (X₂(t) – X₁(t)), где «?» означает знак некоторой операции, а именно — операции умножения. Эта операция дает зависимость одного члена от другого (в нашем случае K от X₂(t) – X₁(t)).

При K = 0, какие бы значения ни принимало выражение X₂(t) – X₁(t), результат K • (X₂(t) – X₁(t)) дает 0. При K = 0.1 значение выражения принудительно ослабляется в 10 раз. То же самое верно и со стороны выражения X₂(t) – X₁(t) — это естественно, ведь общее выражение K • (X₂(t) – X₁(t)) симметрично. Значит, мультипликативная связь моделирует свойство нелинейности и взаимозависимости процессов (один может сводить на нет действие другого). Сложение таким свойством не обладает.

В итоге модель имеет вид:

dX₁(t)/dt = K₂₁ • (X₂(t) – X₁(t)) + K₃₁ • (X₃(t) – X₁(t)).

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

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

dX₂(t)/dt = K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)).

Уравнение изменения температуры атмосферы: dX₃(t)/dt = 0. То есть X₃ = const (X₃ не изменяется).

Уравнение изменения температуры опоры: dX₄(t)/dt = 0. То есть X₄ = const (X₄ не изменяется).

Вся система уравнений в сборе имеет вид:

dX₁(t)/dt = K₂₁ • (X₂(t) – X₁(t)) + K₃₁ • (X₃(t) – X₁(t));

dX₂(t)/dt = K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t));

dX₃(t)/dt = 0;

dX₄(t)/dt = 0.

По физическим соображениям ясно: сколько тепла вытекает из A в B, столько же тепла поступает в В из А, то есть K₂₁ = K₁₂.

К записи общей модели остается добавить конкретные значения коэффициентов теплопроводности (K₁₂ = K₂₁, K₃₁, K₃₂, K₄₂) и начальное состояние системы:

X₁(0) = a,

X₂(0) = b,

X₃(0) = c,

X₄(0) = d,

где a, b, c, d — числа, указывающие температуру соответствующего объекта в момент времени t = 0.

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

СИСТЕМA — потому что в наличии имеется несколько уравнений.

ОБЫКНОВЕННЫХ — производная использована обычная, а не частная, так как использована одна переменная времени.

ДИФФЕРЕНЦИАЛЬНЫХ — в уравнении встречается выражение производной dX/dt.

УРАВНЕНИЙ — в выражениях имеется знак уравнивания.

В КАНОНИЧЕСКОМ ВИДЕ — производная не стоит под знаком какой-либо функции, встречается один раз и только в левой части уравнения. В правой части какие-либо производные отсутствуют.

ЕДИНСТВЕННОЙ — переменная, по которой берут производную, одна во всех уравнениях (время t).

НЕЗАВИСИМОЙ ПЕРЕМЕННОЙ — переменная t не зависит более ни от каких переменных и изменяется сама по себе.

ЛЮБУЮ ДИНАМИЧЕСКУЮ СИСТЕМУ МОЖНО ПРИВЕСТИ К ДАННОМУ ВИДУ!

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

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

dX₁(t)/dt + dX₂(t)/dt + dX₃(t)/dt + dX₄(t)/dt

или

d(X₁(t) + X₂(t) + X₃(t) + X₄(t))/dt

или

dX_{системы}(t)/dt.

В правой части мы имеем следующее:

K₂₁ • (X₂(t) – X₁(t)) + K₃₁ • (X₃(t) – X₁(t)) + K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t))

или (с учетом того, что K₂₁ = K₁₂)

K₃₁ • (X₃(t) – X₁(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)).

И, наконец, вместе:

dX_{системы}(t)/dt = K₃₁ • (X₃(t) – X₁(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)),

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

В случае закрытой системы уравнения имели бы вид:

dX₁(t)/dt = K₂₁ • (X₂(t) – X₁(t)) + K₃₁ • (X₃(t) – X₁(t));

dX₂(t)/dt = K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t));

dX₃(t)/dt = K₁₃ • (X₁(t) – X₃(t)) + K₂₃ • (X₂(t) – X₃(t)) + K₄₃ • (X₄(t) – X₃(t));

dX₄(t)/dt = K₂₄ • (X₂(t) – X₄(t)) + K₃₄ • (X₃(t) – X₄(t)).

При сложении это дает: dX_{системы}(t)/dt = 0 — уравнение закрытой системы, в которую извне нечего не притекает и из которой ничего не истекает. Хотя внутри системы происходят изменения (перераспределение тепла), суммарная температура системы неизменна. Ни первая модель (модель открытой системы), ни вторая модель (модель закрытой системы) не являются «плохими». Просто достигнуты различные цели, даны различные представления исследователя о процессах, о системе в целом.

Динамическая система, которую мы рассматривали выше, — это общий случай. Весьма важно, что из нее всегда можно «даром» получить статическую систему, для чего нужно потребовать: X_i(t + Δt) = X_i(t), в результате чего имеем: dX_i(t)/dt = 0. Другими словами, запись X_i(t + Δt) = X_i(t) означает, что прошлое равно настоящему, то есть состояние системы не меняется. Тогда от уравнений остается:

K₂₁ • (X₂(t) – X₁(t)) + K₃₁ • (X₃(t) – X₁(t)) = 0

K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)) = 0.

Это уравнения статики.

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

Рисунок - Система взаимодействующих тел со слабым нагревателем в задаче теплопроводности

Граф изменится, так как на тело B будет действовать дополнительный объект — нагреватель. Обозначим его температуру переменной X₅(t).

Рисунок - Граф зависимости переменных модели со слабым нагревателем

Изменения в записи будут касаться только второго уравнения, описывающего изменение переменной X₂ во времени. Кроме того, нам понадобится уравнение нагревателя, так как в системе появилась дополнительная переменная X₅(t) (вершина графа) и закон ее изменения должен быть определен.

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

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

Так как такой источник энергии может не только нагревать другие тела, но и нагреваться от них сам, он называется слабым. То есть от услуг такого источника энергии можно отказаться, поскольку он находится в равном с другими телами положении, может влиять на них и может испытывать их влияние. На сильный источник другие тела воздействовать не могут, а сам он — может на них воздействовать, НАВЯЗЫВАЯ ИМ СВОЮ ВОЛЮ БЕЗУСЛОВНО.

Понятно, что в уравнении для тела B появится дополнительное слагаемое, аналогичное ранее рассмотренным:

dX₂(t)/dt = K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)) + K₅₂ • (X₅(t) – X₂(t)).

Уравнение нагревателя может быть записано так:

dX₅(t)/dt = (K • U(t) – X₅(t))/L + K₂₅ • (X₂(t) – X₅(t)),

где U(t) — напряжение питания источника.

Начальные условия X₅(0) и функция работы источника U(t) должны быть заданы:

X₅(0) = h,

U(t) = f(t).

Уравнение связывает причину (напряжение источника питания) и следствие (температура нагревателя). Свойства этого элемента таковы, что при выходе на рабочий режим нагреватель поддерживает постоянную температуру. При включении нагревателя рабочая температура устанавливается со временем, постепенно нарастая, при выключении источника нагреватель постепенно остывает. Очевидно, что процесс нагрева и остывания нагревателя инерционный. Фактически, для описания этих свойств нагревателя достаточно записать апериодический закон (который мы ранее уже обсуждали изменения его температуры X₅(t) по отношению ко входу U(t), что мы и сделали. В уравнении учтен коэффициент усиления K между входом U(t) и выходом X₅(t), инерционность процесса нагрева L и коэффициент теплопередачи K₂₅ между телами X₂(t) и X₅(t).

Теперь рассмотрим гипотезу о сильном нагревателе. Это случай, когда нагреватель устроен таким образом, что излученная им энергия не может отразиться телом B, или если энергия от тела B не может попасть к нагревателю, или если тело B не может «отказаться» от услуг нагревателя.

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

На рисунке изображен граф зависимостей для данного случая.

Рисунок - Граф зависимости переменных модели с сильным нагревателем

Обратите внимание: неважно, нагрет ли нагреватель сильнее тела B или нет — тело B будет получать от него энергию в любом случае, причем, в том объеме, в каком эту энергию будет выделять нагреватель. В уравнении тела B появится дополнительное слагаемое:

dX₂(t)/dt = K₁₂ • (X₁(t) – X₂(t)) + K₃₂ • (X₃(t) – X₂(t)) + K₄₂ • (X₄(t) – X₂(t)) + K₅₂ • X₅(t).

Уравнение нагревателя может быть записано так:

dX₅(t)/dt = f(t).

Начальные условия X₅(0) и функция работы источника f(t) должны быть заданы:

X₅(0) = h.

Но вернемся к расчету движения во времени динамической системы, для которой есть все необходимые данные, есть начальное состояние (X_i(0) = const). Можно по формулам (с использованием метода Эйлера) вычислить скорость ее изменения и новое состояние:

новое состояние := старое состояние + скорость • отрезок времени.

Формально эта запись выглядит так:

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. Простейшим методом численного интегрирования дифференциальных уравнений является метод Эйлера.

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

Рассмотрим практически применение метода Эйлера для расчета процесса изменения температур тел системы. Зададим значения коэффициентов модели: K₁₂ = K₂₁ = 0.2, K₃₁ = 0.1, K₃₂ = 0.05, K₄₂ = 0.1. Зададим начальные условия системы (в момент времени t = 0): X₁(0) = 30° C, X₂(0) = 70° C, X₃(0) = 22° C, X₄(0) = 15° C. Выбираем шаг моделирования Δt равный, например, 0.2 с. Примем конечное значение времени моделирования за T_k = 4 с.

Подставим значения коэффициентов:

X₁(t + Δt) = X₁(t) + [0.2 • (X₂(t) – X₁(t)) + 0.1 • (22 – X₁(t))] • Δt

X₂(t + Δt) = X₂(t) + [0.2 • (X₁(t) – X₂(t)) + 0.05 • (22 – X₂(t)) + 0.1 • (15 – X₂(t))] • Δt

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

Таблица расчета изменения значений переменных системы во времени

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

Рисунок - Расчетные траектории поведения системы тел в задаче теплопроводности

Заметим, что на рисунке показана не вся траектория, а только ее часть. Рассчитайте всю траекторию, подумайте и ответьте на следующие вопросы:

• почему графики в итоге при большом времени рассмотрения стремятся к числу 18.5;
• почему значение переменной X1 сначала увеличивается, а потом падает;
• почему графики имеют переломы производных;
• что надо изменить в условиях задачи, чтобы график X₁ все время убывал?

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

Формальная математическая запись

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

При реализации на математической машине среды Stratum эта запись будет дополнена элементами (функциями) ввода и вывода информации и командами управления (например, «стоп»), а часть элементов будет спрятана (например, при задании начальных условий).

Иногда (когда это удобно) начальные условия могут быть помещены внутрь уравнений. Для этого в следующей записи мы воспользуемся дельта-функцией Дирака.

Напомним, что osc2d(T, X₁) — функция, рисующая на экране точку с координатами T и X₁ для каждого момента времени. Выше данная функция приведена условно: в действительности надо использовать несколько функций для настройки окон, масштаба изображения, цвета, толщины, стиля линии и тому подобного. Подробно реализацию двухмерного осциллографа в среде «Stratum-2000» вы можете посмотреть в тексте имиджа OSCSpace2D.

Приведенный выше вариант формальной математической записи подразумевает, что весь код заключен в рамках одного элемента (имиджа, если использовать термины среды «Stratum-2000»). Будет гораздо более наглядно, если распределить код по отдельным элементам, связав их между собой связями (см. рисунок) — в этом случае структура проекта останется обозримой; кроме того, если потребуется изменить что-то в одном из блоков, то такой подход позволит не менять код в остальных блоках. Еще один плюс такого подхода — отдельные «кирпичики» схемы могут быть применены неоднократно в готовом виде в виде копии имиджа.

Рисунок - Схема проекта в среде Stratum-2000, реализующего полную модель теплопроводности системы тел с элементами интерфейса

Алгоритмическая реализация

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

Рисунок - Блок-схема алгоритма, имитирующего теплообмен между телами

Заметьте, что для реализации нам понадобились следующие инструменты:

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

В теле цикла присутствуют:

• счетчик времени;
• условие выхода из цикла для контроля за временем моделирования;
• блок задания внешних воздействий;
• блок расчета приращений;
• блок расчета нового состояния системы;
• блок вывода текущих состояний.

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

Задание 2

Рассмотрим бак, в котором идет химическая реакция 3A+B-> 2C.

Очевидно, что характеристикой динамики химической реакции является концентрация веществ nA, nB, nC в баке. В начале реакции имеет место большая концентрация nA и nB, и маленькая концентрация nС.

Со временем, по мере того вещества А и В вступают в реакцию концентрации nA, nB уменьшаются, а концентрация nC увеличивается.

Скорость реакции существенно зависит от температуры в баке и от скорости перемешивания жидкостей в баке.

Обратите внимание, если вещества А или В в баке нет, то реакция не идет вовсе и вещество С не образуется. То есть концентрация nA, если она равна нулю, как бы выключает химическую реакцию и, конечно, влияет на скорость ее протекания.

В отличие от первого задания, здесь имеет место нелинейность, взаимозависимость переменных nA и nB.

Так как во времени меняются 3 переменных nA, nB, nC, то ожидается, как минимум, три дифференциальных уравнения с параметрами температуры T и скорости перемешивания V.

Расставьте правильно коэффициенты в уравнениях, указывающие на темп выхода из реакции молекул А и молекул В, они уходят в неравных пропорциях (на одну молекулу В приходится в реакции 3 молекулы А). А также обратите внимание на коэффициент, влияющий на темп образования молекул вещества С на единицу молекулы В.

Задание 3

В задании 2 добавьте в систему три крана, первый – для подачи свежего 100%-го «раствора» вещества А, второй – для вещества В, третий (на дне бака) – для удаления всех веществ в равной доле из зоны реакции. Уровень открытия крана – управляемая величина [литр/сек].

Скорректируйте уравнения системы.

Задание 4

Напишите модель боевых действий – армия синих против армии красных.

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

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

Задание 5

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

Количество травы зависит от наличия солнца и воды.

Количество лис зависит от количества зайцев в лесу и размера популяции лис.

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

Задание 6

Напишите уравнения процесса горения. Имеются горючее в заданном количестве, окислитель (воздух) в неограниченном количестве. Учтите, что во время горения температура горючего не меняется. Горючее вначале процесса подогревают до критической температуры А, по достижении которой оно воспламеняется и во время горения поддерживает температуру А. После выжигания всего горючего температура падает до температуры окружающего воздуха. В результате горения образуется зола, масса которой нарастает во время горения.

Задание 7

Напишите уравнения процесса испарения воды в баке, которую нагревают с температурой нагревателя Т_н. Во время кипения температура воды не меняется, меняется ее масса.