Вычислительная техника в инженерной практике

Размещено на http://www.allbest.ru/

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

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

«Южно-Уральский государственный университет»

(национально-исследовательский университет)

Кафедра «Летательные аппараты»

Курсовая работа

«Вычислительная техника в инженерной практике»

Молявкина Т.П.

Челябинск 2017



Введение

Целью семестровой работы является расчет задач:

1) полета тела в атмосфере под действием силы тяги двигателя,

2) нахождения решения дифференциального уравнения второго порядка, с использованием численных методов, запрограммированных на ЭВМ в программном комплексе Matlab, определение точности полученного решения, устойчивости использованных методов.

Для достижения поставленной цели необходимо:

1) Решить задачу полета тела в атмосфере под действием силы тяги двигателя методом конечных разностей первого, 2-го, 4-го порядков точности, модифицированным методом Эйлера, методом Рунге-Кутты 2-го и 4-го порядков точности, определить точность полученных решений, выявить влияние шага интегрирования на точность получаемого решения.

2) Решить дифференциальное уравнение второго порядка методом прогонки, сравнить результаты с аналитическим решением, выявить влияние шага интегрирования на точность получаемого решения.



Задание №1. Полет материальной точки в атмосфере под действием силы тяги двигателя.

Таблица 1.

№	Исходные данные
	Масса, кг	Площадь, м2	Сила лобового сопротивления	Тяга, Н	(t), град	V(0), м/с	X(0), м	Y(0), м	б(0), град
9	9	0,9	0,9	9		60	80	0,9	40

Дифференциальные уравнения движения тела.

1.1 Метод конечных разностей

1) Дискретный аналог первой производной метода конечных разностей первого порядка точности:



Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:

где

сила лобового сопротивления,

коэффициент лобового сопротивления,

плотность воздуха,

тяга двигателя,

масса тела,

угол наклона траектории (угол атаки).

угол установки двигателя,

б) перемещения вдоль осей:

2) Дискретный аналог первой производной метода конечных разностей второго порядка точности:



Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:

б) перемещения вдоль осей:

3) Дискретный аналог первой производной метода конечных разностей четвертого порядка точности:

Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:



б) перемещения вдоль осей:

1.2 Модифицированный метод Эйлера

Дискретный аналог:

Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:

б) перемещения вдоль осей:



1.3 Метод Рунге-Кутты

Дискретный аналог по методу Рунге-Кутты 4-го порядка:

Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:

б) перемещения вдоль осей:



Дискретный аналог по методу Рунге-Кутты 2-го порядка:

Дискретный аналог для уравнений, описывающих движение тела:

а) проекции скоростей:

б) перемещения вдоль осей:

1.4 Результаты расчета

На рисунке 2 представлены траектории полета тела, полученные численными методами решения ДУ (дифференциальных уравнений) движения тела: МКР (метод конечных разностей), модифицированным методом Эйлера, методом Рунге-Кутты 4-го порядка и методом Рунге-Кутты 2-го порядка, при шаге с.



Рис.2 Траектория полета тела.

1.5 Определение точности полученного решения

Считаем, что метод Рунге-Кутты наиболее точный. Настраиваем его до максимальной точности подбирая число узлов, при котором ошибка меньше 0,001. Метод Рунге-Кутты сходится при шаге 0,002. Но метод МКР и Эйлера при этом значении не сходятся. Поэтому настраиваем шаг.

На рисунках 5 - 8 представлены изменения абсолютных «ошибок» вычислений при шаге с, полученных численными методами решения ДУ (дифференциальных уравнений) движения тела: МКР (метод конечных разностей), модифицированным методом Эйлера, методом Рунге-Кутты 4-го порядка и методом Рунге-Кутты 2-го порядка.



Рис. 3 Абсолютная «ошибка» вычислений МКР.

Рис. 4. Абсолютная «ошибка» вычислений модифицированного метода Эйлера.



Рис. 5 Абсолютная «ошибка» вычислений метода Рунге-Кутты 2-го порядка.

Наиболее точным методом является метод Рунге-Кутты 4-го порядка точности, т.к. при уменьшении шага траектории, полученные при помощи дискретных аналогов методов МКР, модифицированного метода Эйлера, и метода Рунге-Кутты 2-го порядка приближаются к кривой, полученной методом Рунге-Кутты 4-го порядка точности. Абсолютные «ошибки» вычислений всех использованных в расчете методов при уменьшении шага стремятся к нулю.



Задание № 2. Решение ДУ методом «прогонки».

Таблица 2.

Исходные данные
Уравнение	L, м	Y(0)	Y(L)
	9	0,9	9

2.1 Аналитическое решение

Точное решение уравнения имеет вид:

Дискретный аналог для второй производной получим, используя метод конечных разностей первого порядка точности:

тогда

Решение по методу прогонки имеет вид:



и определяются по выражениям:

где

2.2 Результаты расчета

На рисунке 9 представлены графики искомой функции , полученные аналитически и методом «прогонки» с количеством участков разбиения n=100.

Рис. 6. Графики искомой функции , полученные аналитически и методом «прогонки» с количеством участков разбиения n=1000.



2.3 Определение точности полученного решения

Для оценки точности полученного решения определим абсолютную погрешность:

где решение, полученное методом «прогонки»,

аналитическое решение.

полет двигатель программный дифференциальный

Рис. 7. Зависимость абсолютной погрешности решения при изменении количества участков разбиения от n=1000 до n=10000.

Метод прогонки позволяет численно получить решение ДУ второго порядка, при этом точность решения пропорциональна количеству участков разбиения. Абсолютная погрешность решения стремится к нулю при увеличении количества участков разбиения до бесконечности.



Приложение А

clear all, close all

% Исходные данные

dt = 0.000005; % шаг интегрирования, с

m = 9; % масса, кг

S = 0.9; % площадь миделя, м:2

Cx = 0.9; % коэф. лоб. сопро.

t(1) = 0; % начальное время, с

Alpha(1) = 40*pi/180; % угол атаки, град / рад

betta(1)=3*t(1)*pi/180; % угол установки двигателя / рад

x(1) = 80; % нач. координата х, м.

y(1) = 0,9; % нач. координата y, м.

V0 = 60; % нач. скорость, м/с.

Vx(1)= V0 * cos (Alpha(1)); % нач. скорость Vx

Vy(1)= V0 * sin (Alpha(1)); % нач. скорость Vy

YY = 0; % критерий выхода из цикла

i = 0; % счетчик

Ro = 1.205; % плотностьб кг/м3

g = 9.8; % ускорение своб. пад.,м/с2

T= 9; % тяга двигателя,H

% Модифицированный метод Эйлера

xmT(1) = 80;

ymT(1) = 0,9;

x_mT(1) = 80;

y_mT(1) = 0,9;

x_m(1) = 80;

y_m(1) = 0,9;

ax(1)=0;

ay(1)=0;

Vxm(1)=V0 * cos (Alpha(1));

Vym(1)= V0 * sin (Alpha(1));

VxmT(1) = Vxm(1); % тильда

VymT(1)= Vym(1);

AlphaT(1) = 40*pi/180;

bettam(1)=3*t(1)*pi/180;

Alpham(1) = 40*pi/180;

tm(1) = 0;

% Метод Рунге-Кутты 4-го порядка

x_r(1) = 80;

y_r(1) = 0,9;

Alpha1(1) = 40*pi/180;

Alpha2(1) = 40*pi/180;

Alpha3(1) = 40*pi/180;

Vxr(1)=V0 * cos (Alpha(1));

Vyr(1)= V0 * sin (Alpha(1));

X(1) = Cx*Ro*S*(Vxr(1)^2+Vyr(1)^2)/2;

kx0(1)=0;

ky0(1)=0;

kx1(1)=0;

ky1(1)=0;

kx2(1)=0;

ky2(1)=0;

kx3(1)=0;

ky3(1)=0;

bettar(1)=3*t(1)*pi/180;

Alphar(1) = 40*pi/180;

tr(1) = 0;

x_r2(1) = 80;

y_r2(1) = 0,9;

Alphar2(1) = 40*pi/180;

Alphar_2(1) = 40*pi/180;

tr2(1) = 0;

bettar2(1)=3*t(1)*pi/180;

bettar_2(1)=3*t(1)*pi/180;

kx(1)=0;

ky(1)=0;

Vxr2(1)=V0 * cos (Alpha(1));

Vyr2(1)= V0 * sin (Alpha(1));

Error1(1)=0;

SError1(1)=0;

Error2(1)=0;

SError2(1)=0;

%Метод Рунге-Кутты 4го порядка

YY7=y_r2(1);

while (YY7 >=0) && (i<3000000)

%dt=0.005;

i = i + 1;

tr2(i+1) = tr2(i)+dt;

bettar2(i+1)=3*tr2(i+1)*pi/180;

bettar_2(i+1)=3*(tr2(i+1)+dt/2)*pi/180;

Vxr2(i+1) = Vxr2(i)+(-(Cx*Ro*S*((Vxr2(i)+kx(i)/2)^2+(Vyr2(i)+ky(i)/2)^2)/2)*cos(Alphar_2(i))+T*cos(bettar_2(i)))*dt/m;

Vyr2(i+1) = Vyr2(i)+(-(Cx*Ro*S*((Vxr2(i)+kx(i)/2)^2+(Vyr2(i)+ky(i)/2)^2)/2)*sin(Alphar_2(i))-m*g+T*sin(bettar_2(i)))*dt/m;

if Vxr2(i+1)<=0

Alphar2(i+1) = (180*pi/180)+atan(Vyr2(i+1)/Vxr2(i+1));

else

Alphar2(i+1) = atan(Vyr2(i+1)/Vxr2(i+1));

end

kx(i+1)=(-(Cx*Ro*S*((Vxr2(i+1))^2+(Vyr2(i+1))^2)/2)*cos(Alphar2(i+1))+T*cos(bettar2(i+1)))*dt/m;

ky(i+1)=(-(Cx*Ro*S*((Vxr2(i+1))^2+(Vyr2(i+1))^2)/2)*sin(Alphar2(i+1))-m*g+T*sin(bettar2(i+1)))*dt/m;

if (Vxr2(i+1)+kx(i+1))<=0

Alphar_2(i+1) = (180*pi/180)+atan((Vyr2(i+1)+ky(i+1))/(Vxr2(i+1)+kx(i+1)));

else

Alphar_2(i+1) = atan((Vyr2(i+1)+ky(i+1))/(Vxr2(i+1)+kx(i+1)));

end

x_r2(i+1) = x_r2(i) + Vxr2(i)*dt;

y_r2(i+1) = y_r2(i) + Vyr2(i)*dt;

YY7=y_r2(i+1);

end

NN=i;

YYr2 = 0; % критерий выхода из цикла

i = 0; % счетчик

%Метод МКР 1го порядка

m = 9; % масса, кг

S = 0.9; % площадь миделя, м:2

Cx = 0.9; % коэф. лоб. сопро.

t(1) = 0; % начальное время, с

Alpha(1) = 40*pi/180; % угол атаки, град / рад

betta(1)=3*t(1)*pi/180; % угол установки двигателя / рад

x(1) = 80; % нач. координата х, м.

y(1) = 0,9; % нач. координата y, м.

V0 = 60; % нач. скорость, м/с.

Vx(1)= V0 * cos (Alpha(1)); % нач. скорость Vx

Vy(1)= V0 * sin (Alpha(1)); % нач. скорость Vy

YY = 0; % критерий выхода из цикла

i = 0; % счетчик

Ro = 1.205; % плотностьб кг/м3

g = 9.8; % ускорение своб. пад.,м/с2

T= 9; % тяга двигателя,H

%while (YY >=0) && (i<16959)

while (i<NN)

i = i + 1;

t(i+1) = t(i)+dt;

X(i) = Cx*Ro*S*(Vx(i)^2+Vy(i)^2)/2;

Vx(i+1) = Vx(i)+(-X(i)*cos(Alpha(i))+T*cos(betta(i)))*dt/m;

Vy(i+1) = Vy(i)+(-X(i)*sin(Alpha(i))-m*g+T*sin(betta(i)))*dt/m;

x(i+1) = x(i) + Vx(i)*dt;

y(i+1) = y(i) + Vy(i)*dt;

betta(i+1)=3*t(i)*pi/180;

if Vx(i+1)<=0

Alpha(i+1) = (180*pi/180)+atan(Vy(i+1)/Vx(i+1));

else

Alpha(i+1) = atan(Vy(i+1)/Vx(i+1));

end

t(i+1) = t(i) + dt;

YY = y(i+1);

Error1(i+1)= abs(y(i+1)-y_r2(i+1));

end

YYm = 0; % критерий выхода из цикла

i = 0; % счетчик

%Метод Эйлера

%while (YYm >=0) && (i<16959)

while (i<NN)

i = i + 1;

tm(i+1) = tm(i)+dt;

Xm(i) = Cx*Ro*S*(Vxm(i)^2+Vym(i)^2)/2;

ax(i)=-Xm(i)*cos(Alpham(i))+T*cos(bettam(i));

ay(i)=-Xm(i)*sin(Alpham(i))-m*g+T*sin(bettam(i));

VxmT(i+1) = Vxm(i)+ax(i)*dt/m; % тильда

VymT(i+1) = Vym(i)+ay(i)*dt/m; % тильда

bettam(i+1)=3*tm(i)*pi/180;

if VxmT(i+1)<=0

AlphaT(i+1) = (180*pi/180)+atan(VymT(i+1)/VxmT(i+1));

else

AlphaT(i+1) = atan(VymT(i+1)/VxmT(i+1));

end

Vxm(i+1) = Vxm(i)+(ax(i)-(Cx*Ro*S*(VxmT(i)^2+VymT(i)^2)/2)*cos(AlphaT(i+1))+T*cos(bettam(i+1)))*dt/(2*m);

Vym(i+1) = Vym(i)+(ay(i)-(Cx*Ro*S*(VxmT(i)^2+VymT(i)^2)/2)*sin(AlphaT(i+1))-m*g+T*sin(bettam(i+1)))*dt/(2*m);

if Vxm(i+1)<=0

Alpham(i+1) = (180*pi/180)+atan(Vym(i+1)/Vxm(i+1));

else

Alpham(i+1) = atan(Vym(i+1)/Vxm(i+1));

end

x_mT(i+1) = x_m(i) + Vxm(i)*dt; % тильда

y_mT(i+1) = y_m(i) + Vym(i)*dt; % тильда

x_m(i+1) = x_m(i) + (Vxm(i))*dt;

y_m(i+1) = y_m(i) + (Vym(i))*dt;

YYm = y_m(i+1);

Error2(i+1)= abs(y_m(i+1)-y_r2(i+1));

end

YYr = 0; % критерий выхода из цикла

i = 0; % счетчик

%Метод Рунге-Кутты 2-го порядка

%while (YYr >=0) && (i<16959)

while (i<NN)

i = i + 1;

tr(i+1) = tr(i)+dt;

Vxr(i+1) = Vxr(i)+(kx0(i)+2*kx1(i)+2*kx2(i)+kx3(i))/6;

Vyr(i+1) = Vyr(i)+(ky0(i)+2*ky1(i)+2*ky2(i)+ky3(i))/6;

Xr(i+1) = Cx*Ro*S*(Vxr(i+1)^2+Vyr(i+1)^2)/2;

if Vxr(i+1)<=0

Alphar(i+1) = (180*pi/180)+atan(Vyr(i+1)/Vxr(i+1));

else

Alphar(i+1) = atan(Vyr(i+1)/Vxr(i+1));

end

bettar(i+1)=3*(tr(i+1))*pi/180;

kx0(i+1)=dt*(-Xr(i+1)*cos(Alphar(i+1))+T*cos(bettar(i+1)))/m;

ky0(i+1)=dt*(-Xr(i+1)*sin(Alphar(i+1))-m*g+T*sin(bettar(i+1)))/m;

if (Vxr(i+1)+kx0(i+1)/2)<=0

Alpha1(i+1) = (180*pi/180)+atan((Vyr(i+1)+ky0(i+1)/2)/(Vxr(i+1)+kx0(i+1)/2));

else

Alpha1(i+1) = atan((Vyr(i+1)+ky0(i+1)/2)/(Vxr(i+1)+kx0(i+1)/2));

end

kx1(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx0(i+1)/2)^2+(Vyr(i+1)+ky0(i+1)/2)^2)/2)*cos(Alpha1(i+1))+T*cos(3*(tr(i+1)+dt/2)*pi/180))/m;

ky1(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx0(i+1)/2)^2+(Vyr(i+1)+ky0(i+1)/2)^2)/2)*sin(Alpha1(i+1))-m*g+T*sin(3*(tr(i+1)+dt/2)*pi/180))/m;

if (Vxr(i+1)+kx1(i+1)/2)<=0

Alpha2(i+1) = (180*pi/180)+atan((Vyr(i+1)+ky1(i+1)/2)/(Vxr(i+1)+kx1(i+1)/2));

else

Alpha2(i+1) = atan((Vyr(i+1)+ky1(i+1)/2)/(Vxr(i+1)+kx1(i+1)/2));

end

kx2(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx1(i+1)/2)^2+(Vyr(i+1)+ky1(i+1)/2)^2)/2)*cos(Alpha2(i+1))+T*cos(3*(tr(i+1)+dt/2)*pi/180))/m;

ky2(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx1(i+1)/2)^2+(Vyr(i+1)+ky1(i+1)/2)^2)/2)*sin(Alpha2(i+1))-m*g+T*sin(3*(tr(i+1)+dt/2)*pi/180))/m;

if (Vxr(i+1)+kx2(i+1))<=0

Alpha3(i+1) = (180*pi/180)+atan((Vyr(i+1)+ky2(i+1))/(Vxr(i+1)+kx2(i+1)));

else

Alpha3(i+1) = atan((Vyr(i+1)+ky2(i+1))/(Vxr(i+1)+kx2(i+1)));

end

kx3(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx2(i+1))^2+(Vyr(i+1)+ky2(i+1))^2)/2)*cos(Alpha3(i+1))+T*cos(3*(tr(i+1)+dt)*pi/180))/m;

ky3(i+1)=dt*(-(Cx*Ro*S*((Vxr(i+1)+kx2(i+1))^2+(Vyr(i+1)+ky2(i+1))^2)/2)*sin(Alpha3(i+1))-m*g+T*sin(3*(tr(i+1)+dt)*pi/180))/m;

x_r(i+1) = x_r(i) + Vxr(i)*dt;

y_r(i+1) = y_r(i) + Vyr(i)*dt;

YYr = y_r(i+1);

Error3(i+1)= abs(y_r(i+1)-y_r2(i+1));

end

% disp('последняя точка')

% disp([Vxr2(i+1),Vyr2(i+1)])

%Блок вычислений модулей абсолютных погрешностей

%Блок построения графиков

figure;

plot(x,y,'-.',x_r,y_r,'-',x_r2,y_r2,'--',x_m,y_m,'g');

legend('МКР', 'Р-К 2го','Р-К 4го','Метод Эйлера')

%Блок ошибок

figure;

plot(x, Error1);

figure;

plot(x_m, Error2);

figure;

plot(x_r, Error3);

figure(7);

plot(x_r2,y_r2);



Приложение Б

Листинг решения ДУ методом прогонки.

%Задача №2 - решение методом прогонки уравнения d2y/dx2 = x^3

clc

clear all

for k=3:5

N=10^k

%i=1:100;

h=9/(N-1); %Шаг h=yk/(N-1)

y0=0.9;% начальное значение

yk=9;% конечное значение

x=0:h:9; %x=a1:h:b1

a=1;% ничего не меняй

b=-2;% ничего не меняй

c=1;% ничего не меняй

d=@(x)(x^3*h^2); % правая часть уравнения умноженная на h^2

A(1)=(-c)/b;% ничего не меняй

B(1)=(d(x(2))-y0*a)/b;% ничего не меняй

xa=0:h:9 % твой интервал изменения x

y(1)=y0 % ничего не меняй

for i=2:(N-1)

e=a*A(i-1)+b;

A(i)=(-c)/e;

B(i)=(d(x(i+1))-a*B(i-1))/e;

end

y(N)=9;% конечное значение ук

for i=(N-1):-1:2

y(i)=A(i-1)*y(i+1)+B(i-1);

end

%аналитика

for i=1:N

ya(i)=(1/20)*(x(i)^5)-327.15*x(i)+0.9; % уравнение решённое аналитически

yd(i)=abs(ya(i)-y(i));

end

hold on

% subplot(2,1,1)

% plot(x,y,x,ya)

% xlabel('x')

% ylabel('y')

% grid on

% subplot(2,1,2)

plot(x,yd)

xlabel('t')

ylabel('error')

grid on

end

Размещено на Allbest.ru

...