تبلیغات
مهندسی مکانیک - پاسخ سیستم یک درجه ازادی با برانگیزش کلی

پاسخ سیستم یک درجه ازادی با برانگیزش کلی

نویسنده :مسعود شمس
تاریخ:دوشنبه 4 مهر 1390-07:57 ب.ظ

یک نیروی متناوب به یک سیستم یک درجه ازادی وارد می شه که از ما پاسخ سیستم رو می خواد.ابتدا با استفاده از برازش منحنی سری فوریه حاکم بر سیستم رو بدست می اوریم.سپس معادله رو تشکیل داده و با کمک دستور ode45 معادله رو حل می کنیم.


lکه در ان m=0.25 , k=2500 , c=10 می باشد.با ضرب مقدار فشار در سطح مقطع نیرو بدست میاد.پس برای زمانهای 0تا1 یک معادله ضرب در مساحت و در زمان 1تا2 هم یک معادله ضرب در مساحت می نویسیم و اون و تکرار می کنیم.سپس سری فوریه نیرو و محاسبه می کنیم. اسکریپت زیر معادله نیرو را رسم می کند.


% @Masoud Shams
t1=0:0.01:1;
y1=50000*t1*0.0020;
t2=1.01:0.01:2;
y2=50000*(2-t2)*0.0020;
t=0:0.01:6.02;
y=[y1 y2 y1 y2 y1 y2];
plot(t,y)
xlabel('time')
ylabel('force')
axis([0 6.02 0 100])

به ادمه مطلب مراجعه نمایید...


حالا با استفاده از برازش منحنی سری فوریه نیرو رو بدست میاریم.

fun=fit(t',y','fourier5')

fun = 

     General model Fourier5:
     per(x) = a0 + a1*cos(x*w) + b1*sin(x*w) + 
              a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w)+ 
              a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w)
     Coefficients (with 95% confidence bounds):
       a0 =       49.74  (49.68, 49.81)
       a1 =      -40.73  (-40.82, -40.64)
       b1 =      0.7078  (0.5314, 0.8842)
       a2 =    0.003733  (-0.08551, 0.09297)
       b2 =  -0.0001298  (-0.08931, 0.08905)
       a3 =       -4.52  (-4.609, -4.431)
       b3 =      0.2358  (0.1333, 0.3383)
       a4 =    0.003728  (-0.08552, 0.09298)
       b4 =  -0.0002596  (-0.08944, 0.08892)
       a5 =      -1.623  (-1.712, -1.534)
       b5 =      0.1414  (0.04726, 0.2355)
       w =       3.125  (3.124, 3.127)


>> plot(t,y,'r');hold on;plot(fun,'b')

معادله سیستم هم که کاملا مشخصه.معادله مرتبه دو سیستم رو به دو معادله مرتبه اول تبدیل می کنیم.

(m)d2x/dt2+(c)dx/dt+(k)x=f(t)

x(t=0)=0 
dx/dt(t=0)=0 

به عنوان تغییر متغیر داریم:

dx1/dt=x2

و معادله اصلی به صورت مقابل باز نویسی میشود:

mdx2/dt+cx2+kx1=f(t)

حالا معادله دیفرانسیل به صورت دو معادله مرتبه یک در خواهد امد:

dx1/dt=x2
(dx2/dt=1/m(f(t)-cx2-kx1

(f(t همان سری فوریه بدست امده می باشد.حالا نتایج رو در m.file ای ذخیره میکنیم.

function dx=fun(t,x)
dx=zeros(2,1);
a0 =49.74;
a1 =-40.73;
b1 =0.7078;
a2 =0.003733;
b2 =-0.0001298;
a3 =-4.52;
b3 =0.2358;
a4 =0.003728;
b4 =-0.0002596;
a5 =-1.623;
b5 =0.1414;
w =3.125;
m=0.25;
k=2500;
c=10;
dx(1)=x(2);
dx(2)=((1/m)*(a0 + a1*cos(t*w) + b1*sin(t*w) +... 
a2*cos(2*t*w) + b2*sin(2*t*w) + a3*cos(3*t*w) + b3*sin(3*t*w) +... 
a4*cos(4*t*w) + b4*sin(4*t*w) + a5*cos(5*t*w) + b5*sin(5*t*w)-c*x(2)-k*x(1)));

و حالا از تابع ode45 برای حل معادله سیستم استفاده می کنیم.

[t,x]=ode45('fun',[0,6],[0,0]);
figure(1)
plot(t,x(:,1));
hold on
title('Time response of mechanical tranlation system');
xlabel('Time - sec');



خوب پاسخ سیستم رو به صورت یکسری نقاط در اختیار داریم.برای بدست اوردن تابع سیستم دوباره از برازش منحنی استفاده می کنیم:

>> fit(t,x(:,1),'fourier5')
ans = 
     General model Fourier5:
     ans(x) = a0 + a1*cos(x*w) + b1*sin(x*w) + 
              a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w)+ 
              a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w)
     Coefficients (with 95% confidence bounds):
       a0 =     0.01989  (0.01988, 0.01989)
       a1 =    -0.01633  (-0.01634, -0.01632)
       b1 =  8.968e-005  (7.006e-005, 0.0001093)
       a2 = -1.655e-005  (-2.599e-005, -7.109e-006)
       b2 =  8.797e-007  (-9.147e-006, 1.091e-005)
       a3 =   -0.001843  (-0.001852, -0.001833)
       b3 =  3.082e-005  (1.948e-005, 4.216e-005)
       a4 = -1.808e-005  (-2.761e-005, -8.546e-006)
       b4 =  2.016e-006  (-7.793e-006, 1.183e-005)
       a5 =  -0.0006858  (-0.0006954, -0.0006762)
       b5 =  1.991e-005  (9.585e-006, 3.023e-005)
       w =       3.125  (3.124, 3.125)

 >>plot(ans)


اگر شما این مثال رو به روش تحلیلی حل کنید و بکمک متلب اون رو رسم کنید به شکل زیر می رسید.

و اما حل مسئله به روش رانگ کوتا 4 :

a0 =49.74;
a1 =-40.73;
b1 =0.7078;
a2 =0.003733;
b2 =-0.0001298;
a3 =-4.52;
b3 =0.2358;
a4 =0.003728;
b4 =-0.0002596;
a5 =-1.623;
b5 =0.1414;
w =3.125;
m=0.25;
k=2500;
c=10;
dy1=@(t,y1,y2)y2;
dy2=@(t,y1,y2)(1/m)*(a0 + a1*cos(t*w) + b1*sin(t*w) +... 
a2*cos(2*t*w) + b2*sin(2*t*w) + a3*cos(3*t*w) + b3*sin(3*t*w) +... 
a4*cos(4*t*w) + b4*sin(4*t*w) + a5*cos(5*t*w) + b5*sin(5*t*w)-c*y2-k*y1);
h=.01;
t=0:h:10;
y1=zeros(1,length(t));
y2=zeros(1,length(t));
y1(1)=0;
y2(1)=0;
N=length(t);
for i=1:N-1
    k1 = dy1(t(i),y1(i),y2(i));
    m1 = dy2(t(i),y1(i),y2(i));

    k2 = dy1(t(i)+h/2,y1(i)+0.5*k1*h,y2(i)+0.5*m1*h);
    m2 = dy2(t(i)+h/2,y1(i)+0.5*k1*h,y2(i)+0.5*m1*h);

    k3 = dy1(t(i)+h/2,y1(i)+0.5*k1*h,y2(i)+0.5*m2*h);
    m3 = dy2(t(i)+h/2,y1(i)+0.5*k1*h,y2(i)+0.5*m2*h);

    k4 = dy1(t(i)+h,y1(i)*k3*h,y2(i)*m3*h);
    m4 = dy2(t(i)+h,y1(i)*k3*h,y2(i)*m3*h);
    
    y1(i+1) = y1(i) + (h/6)*(k1 + 2*k2 + 2*k3 +k4);
    y2(i+1) = y2(i) + (h/6)*(m1 + 2*m2 + 2*m3 +m4);
end
subplot(211)
plot(t,y1,'r')
grid
ylabel('displacement')
xlabel('time')
subplot(212)
plot(t,y2,'b')
grid
ylabel('velocity')
xlabel('time')



نوع مطلب : اموزش متلب 

داغ کن - کلوب دات کام
نظرات() 
foot complaints
شنبه 18 شهریور 1396 12:01 ب.ظ
I have read so many posts about the blogger lovers but this paragraph is really a fastidious post,
keep it up.
How long does Achilles tendonitis last for?
دوشنبه 16 مرداد 1396 02:57 ب.ظ
Heya! I'm at work browsing your blog from my new
iphone! Just wanted to say I love reading through your blog and look forward to all your posts!
Keep up the excellent work!
How do you get taller in a day?
یکشنبه 15 مرداد 1396 03:13 ب.ظ
Have you ever thought about including a little bit more than just your articles?
I mean, what you say is important and all. However
just imagine if you added some great visuals or videos to give your posts more, "pop"!

Your content is excellent but with pics and videos, this website
could certainly be one of the greatest in its
field. Terrific blog!
manicure
سه شنبه 15 فروردین 1396 07:56 ب.ظ
Hi! Do you use Twitter? I'd like to follow you if that would be okay.

I'm undoubtedly enjoying your blog and look forward to new posts.
 
لبخندناراحتچشمک
نیشخندبغلسوال
قلبخجالتزبان
ماچتعجبعصبانی
عینکشیطانگریه
خندهقهقههخداحافظ
سبزقهرهورا
دستگلتفکر