Методы решения инженерных задач
Задачей моделирования химико-технологических процессов является вычисление характеристик исследуемого объекта в каждой точке аппарата (так называемых полей концентрации, температуры и т. д.). Основное влияние на искомые характеристики оказывает структура потоков - характер движения элементов потоков
в аппарате.
В общем случае для вычисления поля скорости необходимо решить систему дифференциальных уравнений с частными производными второго порядка. Однако записанные уравнения для реальных аппаратов с заданными граничными и начальными условиями в аналитическом виде в большинстве случаев решить не удается. Для решения таких задач используются специализированные комплексы программ вычислительной гидродинамики. Данные программы решают систему дифференциальных уравнений
методами конечных объемов и требуют значительных вычислительных ресурсов.
На практике часто применяются упрощенные модели, в которых изменение характеристик идет только под одной из координат и во времени. Поэтому дифференциальные уравнения в частных производных зачастую встречаются при решении динамических задач, которые зависят от времени.
Лабораторная работа будет связано с изученеим процесса кондуктивного теплообмена. Вывод уравнений кондуктивного теплообмена и теоретическое описание можно посмотреть в пособии
Но алгоритм решения мы рассмотрим на примере задачи из механики: движения струны.
Для решения данной задачи подключим пакеты для решения уравнений:
@time using ModelingToolkit;
@time using MethodOfLines;
@time using OrdinaryDiffEq;
@time using DomainSets;
@time using Plots;
Пакет ModelingToolkit позволяет вводить переменные, и аналитически упрощать систему уравнений перед ее решением. Это может существенно снизить количество уравний и скорость решения.
Пакет MethodOfLines преобразует систему дифференциальных уравнений с частными производными в систему обычных дифференциальных уравнений больешй размерности. Обычно перменной, по которой ищется решение остается время, а градиенты по пространству разбиваются на сетку (одно, двух или трехмерную) и решение ищется в узлах сетки.
Уравнение колебаний струны относится к уравнениям гиперболического типа.
Каждую точку струны можно охарактеризовать значением ее абсциссы . Для определения положения струны в момент времени достаточно знать компоненты вектора смещения точки в момент времени .
Будем предполагать, что смещения струны лежат в одной плоскости (x,u) и что вектор смещения перпендикулярен в любой момент времени к оси x; тогда процесс колебания можно описать одной функцией u(x,t).
Функция u(x,t) характеризует вертикальное перемещение струны. Уравнение колебаний струны:
а=const - зависит от упругости, жесткости, массы и т. д.
Зададим переменные от которых будет зависеть искамая функция:
@parameters x t
Зададим искомую функцию
@variables u(..)
Зададим первую и вторую производную:
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Зададим константу в урванении
a = 5.0
Дифференциальное уравнение или систему уравнений зададим в привычном виде. Знак ~ означает равенство в уравнении.
Можно задать систему уравнений вводя новые уравнения как элементы массива.
eq = [Dxx(u(x,t)) ~ a * Dtt(u(x,t))]
Зададим интервал по времени и пространству (длину струны), на которых будем решать уравнение.
xmin = 0.0
xmax = 6*pi
tmin = 0.0
tmax = 100.0
domains = [x ∈ Interval(xmin, xmax),
t ∈ Interval(tmin, tmax)]
Зададим начальные и граничные условия.
Начальные условия - значения функции U в начальный момент времени. Функция будет завистеть от координаты x. Т.к. мы решаем дифференциальное уравнение второго порядка, то начальных условий должно быть два. Для самой функции u и ее производная - скорость смещения струны
Граничные условия - значения функции, или ее производной на границе диапазона
Допустим в начальный момент времени координата и скорость будет нулевая. И за левый конец струны мы будем ее поднимать и опускать по синусоиде.
bc = [ u(x,0) ~ 0.0, #начальные координаты равны 0
Dt(u(x,0)) ~ 0.0, #начальные скорости равны 0
u(0, t) ~ sin(0.5*t), #условие в начале
]
Создаем проблему с заданной системой уравнений, граничными и начальными условиями
@named pdesys = PDESystem(eq, bc, domains, [x,t], [u(x,t)])
Разделим струну на 100 отрезков
N = 100
discretization = MOLFiniteDifference([x=>N], t)
prob = discretize(pdesys,discretization)
sol = solve(prob, TRBDF2(), saveat=0.1)
Решение будем находить для фиксированных значений и .
discrete_x = sol[x]
solu = sol[u(x, t)]
Построим изменение амплитуды струны по времени.
using Plots
anim = @animate for time = 1:length(discrete_t)
plot(discrete_x, solu[:,time], title="t = $(discrete_t[time])", ylims = (-2.5,2.5), label="") #
xlabel!("x")
ylabel!("u")
end
gif(anim, "struna.gif", fps = 50)
