Моделирование падения аэрозольной частицы в воздухе при наличии сопротивления

Министерство образования и науки Российской Федерации

Новосибирский Государственный Технический Университет

Кафедра ГДУ

РГР по дисциплине «Математическое моделирование физических процессов»

Выполнил: студент гр. Мб-91                                                Проверила: Ивания С.П.

Трушкина М. И.

Новосибирск 2011

Задание 1: Моделирование движения тела в среде при наличии сопротивления.

промоделировать падение аэрозольной частицы в воздухе с высоты Н=100 м, если сила сопротивления воздуха пропорциональна скорости (использовать формулу Стокса), радиус частицы R=50 мкм, плотность материала частицы r = 3,7 10³ кг/м³.

1.Определить скорость падения частицы в момент ее приземления.

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

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

  Алгоритм численного решения:

1.1

 

По первому закону Ньютона: ma=mg- Fсопр

Fсопр=6µπrv

отсюда ускорение будет равно: a=g-6πµrv/m

m=4πρr³/3

µ=1,82 10^-5 кг/м

ρ=3,7 10³ кг/м³

 r = 5 10^-5 м

a= g- 18πµrv/4πr³ρ = g- 9µv/2r²ρ (1)

dv/dt= a где а равно выражению (1).

dy/dt=v

Запишем начальные условия:

V_t=0=0;

Y_t=0=0.

Разобьем область изменения переменной t на N интервалов с шагом Dt и заменим в уравнениях дифференциальные выражения их конечно-разностными приближениями, после чего получим итерационные формулы для скорости в точке сетки:

a_i= g - 9µv/ 2r²ρ;

v_i+1 = v_i + a_idt;

y_i+1 = y_i + v_idt .

Определим предельную скорость при падении частицы, если сила сопротивления определяется формулой Стокса. Сила тяжести равна Vgr, где V, r - объем и плотность материала частицы; подъемная сила равна Vgr^‘ , где r^‘  - плотность среды. Подъемная сила и сила сопротивления направлены противоположно скорости падения, а для установившегося движения сумма всех действующих сил равна нулю. Отсюда

Vg(r - r^‘ ) = 6pmrv_пр

Подставив выражение для объема частицы V = 4/3pr³, получим

v_пред = 2/9r r²g/m.         

Текст программы:

program padenie;

uses crt;

const                     {Запишем в константы величины с которыми мы будем работать}

m=0.0000182;

ro=3700;  

r=0.00005;                       {Радиус}

g=9.8;                   {Ускорение свободного падения}

y0=0;                    {Начальные условия}

v0=0;

h=100;      {Высота}

dt=0.1;      {Изменение времени}

var

a,v,y,t, vpr: real;

f:text;

begin

y:=y0;                   {Присваиваем начальные условия переменным}

v:=v0;

repeat

v:=v+(g-(9*m*v)/(2*r*r*ro))*dt; {Применяем метод Эйлера}

y:=y+v*dt;

t:=t+dt;

vpr:=(2*ro*r*r*g)/(9*m);

Writeln(   'y=', y:0:4,' v=',v:0:4);  {Выводим значения на экран}

writeln('vpr=',vpr:0:4);

until (y>=h);   {Цикл продолжается пока значения координаты y}

readln;       

end.

 Анализ результатов:

Из полученных данных следует что скорость в момент падения частицы равна предельной скорости. v=1,1068 м/с.