БЕЛОРУССКИЙ НАЦИОНАЛЬНЫЙТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ
РЕСПУБЛИКАНСКИЙИНСТИТУТ ИННОВАЦИОННЫХ ТЕХНОЛОГИЙ
КАФЕДРАИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
Курсовая работа
Дисциплина «Математическоемоделирование»
Тема: «Моделированиедвижения парашютиста»
Минск 2008
Содержание
Введение
1. Свободное падение тела с учетом сопротивления среды
2. Формулировка математической модели и ее описание.
3. Описание программы исследования с помощью пакета Simulink
4. Решение задачи программным путем
Список использованных источников
/>/>/>Введение
Формулировка проблемы:
Катапультавыбрасывает манекен человека с высоты 5000 метров. Парашют не раскрывается,манекен падает на землю. Оценить скорость падения в момент удара о землю.Оценить время достижения манекеном предельной скорости. Оценить высоту, накоторой скорость достигла предельного значения. Построить соответствующиеграфики, провести анализ и сделать выводы./>/>/>/>Цель работы:
Научитьсясоставлять математическую модель, решать дифференциальные уравненияпрограммными средствами (используется язык технических вычислений MatLAB 7.0, пакет расширенияSimulink) и анализировать полученные данные о математической модели./>
/>1. Свободноепадение тела с учетом сопротивления среды
При реальных физических движениях тел в газовой или жидкостнойсреде трение накладывает огромный отпечаток на характер движения. Каждый понимает,что предмет, сброшенный с большой высоты (например, парашютист, прыгнувший ссамолета), вовсе не движется равноускоренно, так как по мере набора скоростивозрастает сила сопротивления среды. Даже эту, относительно несложную, задачунельзя решить средствами “школьной” физики: таких задач, представляющихпрактический интерес, очень много. Прежде чем приступать к обсуждениюсоответствующих моделей, вспомним, что известно о силе сопротивления.
Закономерности, обсуждаемые ниже, носят эмпирический характер иотнюдь не имеют столь строгой и четкой формулировки, как второй закон Ньютона.О силе сопротивления среды движущемуся телу известно, что она, вообще говоря,растет с ростом скорости (хотя это утверждение не является абсолютным). Приотносительно малых скоростях величина силы сопротивления пропорциональнаскорости и имеет место соотношение, /> где /> определяется свойствамисреды и формой тела. Например, для шарика /> — это формула Стокса, где /> —динамическая вязкость среды, r — радиус шарика. Так, для воздуха при t = 20°С и давлении 1 атм /> = 0,0182 H.c.м-2 для воды 1,002 H.c.м-2, для глицерина 1480H.c.м-2.
Оценим, при какой скорости для падающего вертикально шара силасопротивления сравняется с силой тяжести (в движение станет равномерным).
Имеем
/>
или
/> (1)
Пусть r= 0,1 м, />= 0,8 кг/м (дерево). При падении ввоздухе /> м/с,в воде />17м/с, в глицерине />0,012 м/с.
На самом деле первые два результата совершенно не соответствуютдействительности. Дело в том, что уже при гораздо меньших скоростях силасопротивления становится пропорциональной квадрату скорости: />. Разумеется, линейнаяпо скорости часть силы сопротивления формально также сохранится, но если />, то вкладом /> можнопренебречь (это конкретный пример ранжирования факторов). О величине k2 известно следующее: онапропорциональна площади сечения тела S, поперечного по отношению к потоку, и плотностисреды /> изависит от формы тела. Обычно представляют k2 = 0,5сS/>, где с — коэффициент лобовогосопротивления — безразмерен. Некоторые значения с (для не очень большихскоростей) приведены на рис.1.
Придостижении достаточно большой скорости, когда образующиеся за обтекаемым теломвихри газа или жидкости начинают интенсивно отрываться от тела, значение с внесколько раз уменьшается. Для шара оно становится приблизительно равным 0,1.Подробности можно найти в специальной литературе.
Вернемся куказанной выше оценке, исходя из квадратичной зависимости силы сопротивления отскорости.
Имеем
/>
или
/> (2)
для шарика
/> (3)
/>
Диск
Полусфера
Полусфера
Шар
Каплевидное тело
Каплевидное тело
с = 1,11
с = 1,33
с = 0,55
с = 0,4
с = 0,045
с = 0,01
Рис1. Значения коэффициента лобового сопротивления для некоторых тел,поперечное сечение которых имеет указанную на рисунке форму
Примем r = 0,1 м, />=0,8.103 кг/м3 (дерево). Тогда длядвижения в воздухе (/>= 1,29 кг/м3 ) получаем />/>18 м/с, в воде(/>= 1.103 кг/м3 )/>/>0,65 м/с, вглицерине (/>=1,26.103 кг/м3 ) />/> 0,58 м/с.
Сравнивая с приведенными выше оценками линейной части силысопротивления, видим, что для движения в воздухе и в воде ее квадратичная частьсделает движение равномерным задолго до того, как это могла бы сделать линейнаячасть, а для очень вязкого глицерина справедливо обратное утверждение.Рассмотрим свободное падение с учетом сопротивления среды. Математическаямодель движения — уравнение второго закона Ньютона с учетом двух сил,действующих на тело: силы тяжести и силы сопротивления среды:
/> (4)
Движениеявляется одномерным; проецируя векторное уравнение на ось, направленнуювертикально вниз, получаем
/> (5)
Вопрос, который мы будем обсуждать на первом этапе, таков: каковхарактер изменения скорости со временем, если все параметры, входящие вуравнение (7) заданы? При такой постановке модель носит сугубо дескриптивныйхарактер. Из соображений здравого смысла ясно, что при наличии сопротивления,растущего со скоростью, в какой-то момент сила сопротивления сравняется с силойтяжести, после чего скорость больше возрастать не будет. Начиная с этогомомента, />,и соответствующую установившуюся скорость />можно найти из условия />=0, решая недифференциальное, а квадратное уравнение. Имеем
/> (6)
(второй — отрицательный — корень, естественно, отбрасываем). Итак,характер движения качественно таков: скорость при падении возрастает от /> до />. Как и покакому закону – это можно узнать, лишь решив дифференциальное уравнение (7).
Однако даже в столь простой задаче мы пришли к дифференциальномууравнению, которое не относится ни к одному из стандартных типов, выделяемых вучебниках по дифференциальным уравнениям, допускающих очевидным образоманалитическое решение. И хотя это не доказывает невозможность егоаналитического решения путем хитроумных подстановок, но они не очевидны. Допустим,однако, что нам удастся найти такое решение, выраженное через суперпозициюнескольких алгебраических и трансцендентных функций – а как найти законизменения во времени перемещения? Формальный ответ прост:
/> (7)
но шансы на реализацию этой квадратуры уже совсем невелики. Дело втом, что класс привычных нам элементарных функций очень узок, и совершеннообычна ситуация, когда интеграл от суперпозиции элементарных функций не можетбыть выражен через элементарные функции в принципе. Математики давно расширилимножество функций, с которыми можно работать почти так же просто, как сэлементарными (т. е. находить значения, различные асимптотики, строить графики,дифференцировать, интегрировать). Тем, кто знаком с функциями Бесселя,Лежандра, интегральными функциями и еще двумя десятками других, так называемыхспециальных функций, легче находить аналитические решения задач моделирования,опирающихся на аппарат дифференциальных уравнений. Однако даже получениерезультата в виде формулы не снимает проблемы представления его в виде,максимально доступном для понимания, чувственного восприятия, ибо мало ктоможет, имея формулу, в которой сопряжены логарифмы, степени, корни, синусы итем более специальные функции, детально представить себе описываемый ею процесс- а именно это есть цель моделирования.
В достижении этой цели компьютер — незаменимый помощник.Независимо от того, какой будет процедура получения решения — аналитической иличисленной, — задумаемся об удобных способах представления результатов.Разумеется, колонки чисел, которых проще всего добиться от компьютера (что притабулировании формулы, найденной аналитически, что в результате численногорешения дифференциального уравнения), необходимы; следует лишь решить, в какойформе и размерах они удобны для восприятия. Слишком много чисел в колонке бытьне должно, их трудно будет воспринимать, поэтому шаг, с которым заполняется таблица,вообще говоря, гораздо больше шага, с которым решается дифференциальноеуравнение в случае численного интегрирования, т.е. далеко не все значения /> и />, найденныекомпьютером, следует записывать в результирующую таблицу (табл. 2).
Таблица 2
Зависимость перемещения и скорости падения от времени (от 0 до 15с)t(c) S(m)
/>(м/с) t(c) S(m)
/>(м/с)
1
2
3
4
5
6
7
4.8
18.7
40.1
66.9
97.4
130.3
164.7
9,6
17,9
24,4
28,9
31,9
33,8
35,0
8
9
10
11
12
13
14
15
200.1
235.9
272.1
308.5
345.0
381.5
418.1
454.7
35.6
36.0
36.3
36.4
36.5
36.6
36.6
36.6
Кроме таблицы необходимы графики зависимостей /> и />; по ним хорошо видно,как меняются со временем скорость и перемещение, т.е. приходит качественноепонимание процесса.
Еще один элемент наглядности может внести изображение падающеготела через равные промежутки времени. Ясно, что при стабилизации скоростирасстояния между изображениями станут равными. Можно прибегнуть и к цветовойраскраске — приему научной графики, описанному выше.
Наконец, можно запрограммировать звуковые сигналы, которыеподаются через каждый фиксированный отрезок пути, пройденный телом — скажем,через каждый метр или каждые 100 метров — смотря по конкретным обстоятельствам.Надо выбрать интервал так, чтобы вначале сигналы были редкими, а потом, сростом скорости, сигнал слышался все чаще, пока промежутки не сравняются. Такимобразом, восприятию помогают элементы мультимедиа. Поле для фантазии здесьвелико.
Приведем конкретный пример решения задачи о свободно падающемтеле. Герой знаменитого фильма “Небесный тихоход” майор Булочкин, упав с высоты6000 м в реку без парашюта, не только остался жив, но даже смог снова летать.Попробуем понять, возможно, ли такое на самом деле или же подобное случаетсятолько в кино. Учитывая сказанное выше о математическом характере задачи,выберем путь численного моделирования. Итак, математическая модель выражаетсясистемой дифференциальных уравнений.
/> (8)
Разумеется, это не только абстрактное выражение обсуждаемойфизической ситуации, но и сильно идеализированное, т.е. ранжирование факторовперед построением математической модели произведено. Обсудим, нельзя липроизвести дополнительное ранжирование уже в рамках самой математической моделис учетом конкретно решаемой задачи, а именно — будет ли влиять на полетпарашютиста линейная часть силы сопротивления и стоит ли ее учитывать примоделировании.
Так как постановка задачи должна быть конкретной, мы примемсоглашение, каким образом падает человек. Он опытный летчик и навернякасовершал раньше прыжки с парашютом, поэтому, стремясь уменьшить скорость, онпадает не “солдатиком”, а лицом вниз, “лежа”, раскинув руки в стороны. Ростчеловека возьмем средний — 1,7 м, а полуобхват грудной клетки выберем вкачестве характерного расстояния — это приблизительно 0,4 м. для оценки порядкавеличины линейной составляющей силы сопротивления воспользуемся формулойСтокса. Для оценки квадратичной составляющей силы сопротивления мы должныопределиться со значениями коэффициента лобового сопротивления и площадью тела.Выберем в качестве коэффициента число с=1,2 как среднее между коэффициентамидля диска и для полусферы (выбор дня качественной оценки правдоподобен). Оценимплощадь: S= 1,7 ∙ 0,4 = 0,7(м2).
В физических задачах надвижение фундаментальную роль играет второй закон Ньютона. Он гласит, чтоускорение, с которым движется тело, прямо пропорционально действующей на негосиле (если их несколько, то равнодействующей, т.е. векторной сумме сил) и обратнопропорционально его массе:
/>.
Так для свободнопадающего тела под действием только собственной массы закон Ньютона примет вид:
/>
Или в дифференциальномвиде:
/>
/>
Взяв интеграл от этоговыражения, получим зависимость скорости от времени:
/>
Если в начальный момент V0 = 0, тогда />.
Далееопределим зависимость высоты от времени, для чего проинтегрируем последнеевыражение.
/>
/>
/>.
Выясним, прикакой скорости сравняются линейная и квадратичная составляющие силысопротивления. Обозначим эту скорость /> Тогда
/>
или
/>
Ясно, что практически с самого начала скорость падения майораБулочкина гораздо больше, и поэтому линейной составляющей силы сопротивленияможно пренебречь, оставив лишь квадратичную составляющую.
После оценки всех параметров можно приступить к численному решениюзадачи. При этом следует воспользоваться любым из известных методовинтегрирования систем обыкновенных дифференциальных уравнений: методом Эйлера,одним из методов группы Рунге — Кутта или одним из многочисленных неявныхметодов. Разумеется, у них разная устойчивость, эффективность и т.д. — этисугубо математические проблемы здесь не обсуждаются.
Вычисления производятся до тех пор, пока не опустится на воду.Примерно через 15 с после начала полета скорость становится постоянной иостается такой до приземления. Отметим, что в рассматриваемой ситуациисопротивление воздуха радикально меняет характер движения. При отказе от егоучета график скорости, изображенный на рисунке 2, заменился бы касательной кнему в начале координат.
/>
Рис. 2. График зависимости скорости падения от времени/>2. Формулировкаматематической модели и ее описание
парашютист падениесопротивление математическая модель
Припостроении математической модели необходимо соблюдение следующих условий:
— манекенмассой 50 кг соответственно падают в воздухе с плотностью 1,225 кг/м3;
— на движениевлияют только силы линейного и квадратичного сопротивления;
— площадьсечения тела S=0.4м2;
Тогда длясвободно падающего тела под действием сил сопротивления закон Ньютона приметвид:
/>,
где a – ускорение тела, м/с2,
m – его масса, кг,
g – ускорение свободногопадения на земле, g = 9,8 м/с2,
v – скорость тела, м/c,
k1 – линейный коэффициентпропорциональности, примем k1 = β = 6πμl (μ – динамическаявязкость среды, для воздуха μ = 0,0182 Н.с.м-2; l – эффективная длина,примем для среднестатистического человека при росте 1,7 м и соответствующемобхвате грудной клетки l = 0,4 м),
k2 – квадратичныйкоэффициент пропорциональности. K2 = α = С2ρS. В данном случаедостоверно можно узнать лишь плотность воздуха, а площадь манекена S и коэффициент лобовогосопротивления С2 для него определить сложно, можно воспользоваться полученнымиэкспериментальными данными и принять K2 = α = 0,2.
Тогда получим законНьютона в дифференциальном виде:
/>
Так как
/>
Тогда можно составитьсистему дифференциальных уравнений:
/>
Математическая модель припадении тела в гравитационном поле с учетом сопротивления воздуха выражаетсясистемой из двух дифференциальных уравнений первого порядка./>3. Описаниепрограммы исследования с помощью пакета Simulink
Для имитационного моделирования движения парашютиста в системеMATLAB используем элементы пакета расширения Simulink. Для задания величинначальной высоты — H_n,конечной высоты — H_ k, числа — pi, μ – динамическая вязкость среды — my,обхват — R, массе манекена m, коэффициент лобового сопротивления — c, плотность воздуха — ro, площадь сечения тела — S, ускорение свободногопадения — g,начальная скорость — V_n используем элемент Constant находящийся в Simulink/Sources (рисунок 3).
/>
Рисунок 3. Элемент Constant
Для операцииумножения используем блок Product, находящийся в Simulink/Math Operations/Product (рисунок 4).
/>
Рисунок. 4
Для ввода k1 – линейногокоэффициента пропорциональности и k2 – квадратичного коэффициента пропорциональностииспользуем элемент Gain, находящийся в Simulink/Math Operations/Gain (Рисунок. 5.)
/>
Рисунок. 5
Дляинтегрирования – элемент Integrator. Находящийся в Simulink/Continuous/Integrator. Рисунок. 6.
/>
Рисунок. 6
Для выводаинформации используем элементы Display и Scope. Находящиеся в Simulink/Sinks. (Рисунок. 7)
/>
Рисунок. 7
Математическаямодель для исследования с использованием вышеперечисленных элементов,описывающая последовательный колебательный контур приведена на рисунке 8.
/>
Рисунок. 8
Программа исследований
1. Исследование графиказависимости высоты от времени и скорости от времени масса парашютиста равна50кг.
/>
Рисунок 9
Из графиковвидно, что при расчете падения парашютиста массой 50 кг, следующие данные:максимальная скорость равна 41,6 м/с и время равно 18с, и должна достигатьсячерез 800 м падения, т.е. в нашем случае на высоте около 4200 м.
/>
Рисунок. 10
2.Исследование графика зависимости высоты от времени и скорости от времени массапарашютиста равна 100кг.
/>
Рисунок 11
/>
Рисунок 12
С массойпарашютиста 100 кг.: максимальная скорость равна 58 м/с и время равно 15с, идолжна достигаться через 500 м падения, т.е. в нашем случае на высоте около4500 м. (рисунок.11., рисунок. 12).
Выводы пополученным данным, которые справедливы для манекенов, отличающихся толькомассой, но с одинаковыми размерами, формой, типом поверхности и другимипараметрами, определяющими внешний вид объекта.
Легкийманекен при свободном падении в гравитационном поле с учетом сопротивлениясреды достигает меньшей предельной скорости, но за меньший промежуток времении, естественно, при одинаковой начальной высоте – в более низкой точкетраектории, чем тяжелый манекен.
Чем тяжелееманекен, тем быстрее он достигнет земли./> 4. Решениезадачи программным путем
М-файлфункции parashut.m:
%Функция моделированиядвижения парашютиста
functiondhdt=parashut(t,h)
global k1 k2 g m
% система ДУ первого порядка
dhdt(1,1)= -h(2);
dhdt(2,1)=(m*g-k1*h(2)-k2*h(2)*h(2))/m
М-файл вывода результатов parashutist.m:
% Моделирование движенияпарашютиста
% Васильцов С. В.
clc
global h0 g mk1 k2 a
% k1-линейный коэффициентпропорциональности, определяющийся свойствами среды и формой тела. ФормулаСтокса.
k1=6*0.0182*0.4;
%k2-квадратичныйкоэффициент пропорциональности, пропорционален площади сечения тела,поперечного по
%отношения к потоку,плотности среды и зависит от формы тела.
k2=0.5*1.2*0.4*1.225
g=9.81; % ускорениесвободного падения
m=50; % масса манекена
h0=5000; % высота
[t h]= ode45(@parashut,[0200],[h0 0] )
r=find(h(:,1)>=0);
s=length(r);
b=length(t);
h(s+1:b,:)=[];
t(s+1:b,:)=[];
a=g-(k1*-h(:,2)+k2*h(:,2).*h(:,2))/m% вычисляем ускорение
% Построение графиказависимости высоты от времени
subplot(3,1,1),plot(t,h(:,1),'LineWidth',1,'Color','r'),grid on;
xlabel('t,c'); ylabel('h(t), m');
title('График зависимости высоты от времени', 'FontName','Arial','Color','r','FontWeight','bold');
legend('m=50 kg')
% Построение графиказависимости скорости от времени
subplot(3,1,2),plot(t,h(:,2),'LineWidth',1,'Color','b'),grid on;
xlabel('t,c');
ylabel('V(t),m/c');
Title('График зависимости скорости от времени', 'FontName','Arial','Color','b','FontWeight','bold');
legend('m=50 kg')
% Построение графиказависимости ускорения от времени
subplot(3,1,3),plot(t,a,'-','LineWidth',1,'Color','g'),grid on;
text (145,0,'t, c');
ylabel('a(t),m/c^2');
Title('График зависимости ускорения от времени', 'FontName','Arial','Color','g','FontWeight','bold');
legend('m=50 kg')
Экранная форма вывода графиков.
/>/>/>/>/>/>/>
/>Список использованных источников
1. Вся физика. Е.Н.Изергина. – М.: ООО «Издательство «Олимп», 2001. – 496 с.
2. Касаткин И. Л. Репетиторпо физике. Механика. Молекулярная физика. Термодинамика/ Под ред. Т. В. Шкиль.– Ростов Н/Д: изд-во «Феникс», 2000. – 896 с.
3. Компакт-диск «СамоучительMathLAB». ООО «Мультисофт», Россия, 2005.
4. Методические указания кКурсовой работе: дисциплина Математическое моделирование. Движение тела приучете сопротивления среды. – Минск. РИИТ БНТУ. Кафедра ИТ, 2007. – 4 с.
5. Решение системдифференциальных уравнений в Matlab. Дубанов А.А. [Электронный ресурс]. – Режим доступа: rrc.dgu.ru/res/exponenta/ educat/systemat/dubanov/index.asp.htm;
6. Энциклопедия д.д. Физика.Т. 16. Ч.1. с. 394 – 396. Сопротивление движению и силы трения. А. Гордеев./Глав. ред. В.А. Володин. – М. Аванта+, 2000. – 448 с.
7. Matlab Function Reference [Электронный ресурс]. –Режим доступа: matlab.nsu.ru/Library/Books/Math/MATLAB/help/techdoc/ref/.