Mechatronics Design Project-차 례-1. 서론I. 목 표2. 본론I. 가정 및 제한II. 회로도III. 사용된 부품IV. 회로 분석A. 입력전압 Offset의 변화에 따른 회로의 전압 graphB. Steady state에서의 값 분석C. Capacitor에서 T-V, T-I 분석D. 충전 표시 LED 와 과충전 방지 Zener diode의 역할3. 결론1. 서론I. 목 표- 아래의 조건을 만족시키는 Constant Current Battery Charger 를 PSPICE 를 이용해 디자인 하라.A. 9V Ni-Cd 배터리는 9.6V의 실 공칭전압을 가지며 최소전류 20mA 최대 전류 40mA가 요구된다. 더욱이 완충 되엇을때 배터리의 전압은 10.4V이다.B. 입력 직류 전원에 약간의 Fluction 을 고려해야한다. 이러한 Fluction을 견딜수 있어야 한다.C. BJT 를 사용한 회로를 제작하여야 한다.D. 전류를 정류하기 위해 제너다이오드 혹은 다른 전기 부품을 사용해라E. AC/DC 컨버터는 사용할 필요가 없다. 공급되는 전류는 약간의 Fluction 을 포함한 직류이다.F. 배터리 모델이 필요하다면 큰 Capacitor(s) 와 작은 Resistior(s)를 이용해 만들어서 사용하라.2. 본론I. 가정, 제한- 문제의 전지를 20Ω 의 Resistor 과 100F 의 Capacitor 로 이루어진 공칭전압 9.6V의 단일 전지로 가정하였다. 충전전압이 다른 전지의 직렬연결 충전은 전지의 과충전을 불러일으킬 수 있는데 이 회로에선 그러한 문제에 대해선 고려하지 않았다.- 충전을 요하는 전지가 완전히 방전되었을 경우로 가정하고 충전기를 설계 하였다. 만약 충분히 방전되지 않은 전지를 회로에 충전시킬 경우 충전시간은 단축되어 질수 있지만 Ni-Cd 전지 특성상 0.04~0.08v의 전압강하가 일어날 수 있다.- 실제적인 Ni-Cd 전지의 경우에는 실제 용량의 1/10 비율로 14~16시간동안 충전해 완충 할 수 있지만 이 충전기의 경우 10시간동안 충전할 경우 완충될 수 있게 설계하였다.- 직류전원에서 Capacitor가 포함되어 있을 때, 시간에 따른 전류, 전압의 변화를 Pspice에서 나타낼 수 없기 때문에 (직류 : Capacitor를 Open circuit 으로 간주) Capacitor에 전하가 모두 충전된 후 Steady state에서의 BJT 의 Emitter 전압이 충전하는 동안 일정하게 작용했다고 가정하였다. 이는 Capacitor가 완충되어 회로내에서 Open circuit 취급 당할때의 Emiiter 전압과 같다.- Random Fluction 을 가진 DC 전압을 만들기에 제한이 있어 입력전압은 Offset을 가진 교류전압을 이용하였다. 이때 Voffset > 10Vampl 으로 하여 small fluction 의 크기에 제한을 두었다.- Ni-Cd 전지의 경우엔 -20℃~60℃의 범위에서 사용 가능하여 과열방지에 주의하여야 하지만 이 회로의 경우 과열에 의한 문제에 관하여선 고려하지 않았다.II. 회로도전지부분III. 사용된 부품명칭Part name설 명수 량(EA)R~ResistorIEC 표준계열 공칭저항 사용 R6의 E24계열을 제외하곤 E96계열을 사용하되 정확하지 않은 값의 저항은 가장 비슷한 값으로 찾아서 사용한다.1KΩ x 4160Ω x 1500Ω x 120Ω x 1D1Zener diode9.1Version에서 제공하는 1N750다이오드를 사용하였다. 이는 프로그램 자체에서 디자인되어 제공하는 Zener diode 가 1N750뿐이었기 때문에다. 회로에서 필요한 전압을 만들어 내기 위해서 Vz(Break voltage) 는 30V로 수정하였다. 나머지 특성은 1N750 과 같다.1D7동일한 1N750 다이오드를 사용하면 Vz를 수정할 수 가 없어서 임의로 SKKU 다이오드를 만들어서 사용하였다 Vz = 12.3V 이며 나머지 특성은 1N750과 같다.1D3LEDLED 다이오드를 대체하기 위해 일반 다이오드(1N4002)를 사용하였다 LED 다이오드와 비슷한 전압강하를 내도록 Vj = 1.5V 로 수정하였다. 나머지 특성은 1N4002 와 같다각각 1D5LED 다이오드를 대체하기 위해 일반 다이오드(1N4148)를 사용하였다. D3보다 민감하게 전압에 반응해야하여 Vj = 0.7V 로 수정하여 하용하였다. 나머지 특성은 1N4148 과 같다.Q1BJT2N3904 Transister이다. Vbe = 0.75V이다. 별다른 특성 수정 없이 사용하였다.1C1Capacitor단순한 100F의 Capacitor이다.1IV. 회로 분석A. 입력전압 Offset의 변화에 따른 회로의 전압 graph입력전압 회로도에선 임의로 VOffset을 100V로 표시했지만 회로에 적절한 전압을 찾기 위해서 Voffset의 값을 0V에서 300V 까지 Sweep 시켜 D1 을 거친 이후 안정적인 전류가 흐르는 최소 전압을 찾아야 한다. 이를 확인하기 위한 그래프는 다음과 같다.그래프를 확인해보면 Voffset이 90V 이상일때 Zener diode가 제대로 작동하여 30V의 일정한 전압이 측정됨을 알 수 있다. 그렇기 때문에 회로에 제공될 직류전원은 90V 이상의 전압이여야만 한다. 이때는 Zener diode가 제대로 작동하여 10% 이내의 Fluction을 무시할수 있을정도로 작게 해준다. 90V일때 시간에 따른 전압그래프를 확인해보면 입력전원에선 크게 나타난 진폭이 거의 없어졌음을 확인 할 수 있다.B. Steady state에서의 값 분석앞서 가정한 것 처럼 전지가 있는 회로에 연결된가 Capacitor의 충전 여부와 상관 없이 Steady state의 전압으로 일정하므로 Capacitor에 걸리는 전압을 구하기 위해선 일단결정해야 한다. 이는 BJT 인 Q1과 주변 저항, 입력전압(30V)를 통해 결정할 수 있다. 회로를 보면 D7 에 의해 정류된 30V 전압이 양단에 걸리게 되고 이는 곧 아래 그림과 같은 형태가 된다.이는이 양변에 작용한 회로로 볼수 있으며 좌측을 간략화 하기 위해 Thevinin equivalent 회로로 변환하면 다음과 같다.이때 Rth 값은 137.9Ω Vth 값은 25.9V 이다. 이를 트랜지스터의 β값을 이용해 풀어내려 했으나 다른 회로를 통해 β값을 측정하였으나 전류에 따라 다르게 나오는 걸 확인하고, 그냥 Pspice 시뮬레이션을 통해서를 구하였다. Thevinin eqiovalent 회로와 변환전 회로 둘 모두에서값이 22.8V 가 측정되었으므로 변환을 통해 회로를 분석하려 한 시도가 틀리진 않았다는것 알 수 있다. 위 회로를 Pspice로 구현하면값이 22.8V 가 나오며 회로 확인용으로 집어넣은 LED의 전압강하량 1.5~2V를 제외하면 필요한와 근접한 전압(20.8V±0.5V 수준)이 나오게 된다. 이때의 전압분포를 Pspice로 나타낸 그림은 다음과 같다.이때 D7의 경우 D7의 제너전압보다 약간 큰값을 정류하는데 이는 D7이 이상적인 제너 다이오드가 아니라 정확한 Vz값에서 수직으로 Break 되지 않기 때문으로 생각된다.C. Capacitor에서 T-V, T-I 분석Pspice 로는 Capacitor를 Open circuit 으로 보기 때문에 Pspice로는 시간에 따른 전압 충전 혹은 전류 흐름을 나타낼 수 없다. 이를 분석하기 위해 위의 가정 “Emitter Voltage()는 Steady state에서의 전압으로 Capacitor의 충전 여부와 관계없이 일정하다.”를 사용하여 미분방정식을 풀었다. Capacitor에 전압이 가해졌을때 시간에 따른 Capacitor의 전압와의 관계는 다음과 같다.이때 R7,R6,C1 에 걸리는 전류는로 일정하므로이고, 이를 시간에 관해 미분하여 I 를 t 에 관한 식으로 풀면 다음과 같다.이를에 대입하면이렇게 구해진,에 각각 상수()를 집어넣어 계산하면 t=0일때 최대전류 40mA t=36000 일때 최소전류 20mA 완전충전 전압 10.4V 임을 확인할 수 있다. 이를 그린 그래프는 다음과 같다.위 그래프를 보면 t-I 그래프의 경우에 찍힌 두점은 (0,4),(3.6,2) 인데 이는 각각 t=0s 일때 I(t)=40mA 임을 t=36000s 일때 I(t)=20mA 임을 나타내고 있으며 t-V 그래프의 경우에 표시된 두 직선은 각각 y=10.4 x=3.6 인데 이는 t=36000s 에서 V가 10.4V 임을 의미한다.위 그래프의 분석값을 살펴보면 전지의 충전이 t=0s, I=40mA, V=0V에서 시작해서 t=36000s 일 때 I=20mA, V=10.4V로 완충됨을 알 수 있다. 그래프의 기울기가 너무 급격하게 보이기도 하는데 이는 시간축이 0s-36000s 이 될 경우 구분짓기 어려워 그래프를 축소하는 과정에서 일어난 왜곡 때문에 이렇게 보이는 것이지, 실지로 시간에따른 전압변화는 굉장히 천천히 일어나게 된다.D. 충전 표시 LED 와 과충전 방지 Zener diode의 역할회로도에서 전지 옆을 보면 Zener diode 와 LED 2개가 있음을 확인할 수 있다. 이는 전지의 충전상태를 표시해주는 LED이자 전지의 과충전을 방지해주는 역할을 한다. B 의 분석에서 보면 전지가 충전이 다 되는 때는 Vc(t)가 10.4V일때 이다. 이때의 전류 I(t)는 20mA 이다. 그렇기 때문에 총 건전지에 걸리는 전압은 10.8V가 되게 된다. 그리고 여기에 충전을 표시하기 위한 D3의 전압강하(Vj)가 1.5V 이기 때문에 총 D7과 병렬연결된 회로의 전압은 12.3V가 되게 된다. 완충된 캐패시터보다 Ve가 크기 때문에 만약 D7이 없다면 과충전되어 위험을 초래할 수 있으나 D7이 12.3V가 넘어가게 되면 회로를 D5쪽으로 바꾸게 되고 LED에 불이들어오게 되어 완충되었음을 알림과 동시에 12.3V이상의 전류를 억제하여 전지의 과충전을 막게된다. 충전이 되기 이전 즉 전지에 12.3V 이하의 전류가 가해질 때 에는 D7이 아무런 영향을 주지 않으므로 D3 LED에 불이 들어와 충전이 되고 있음을 알 수 있다.
Fast response 를 만족한는 Ki = 70 값에선 kp와 Kd 가 높아질 수록 벨트 장력값이 심하게 불안정한 모습을 확인 할 수 있었다. 이는 Small oscilliation 을 만족하는 Ki=20 에서의응답에서도 동일한 모습 이었다. 그렇기 때문에 각각의 용도에 따라 빠른 반응을 원한다면Ki 값을 적당히 높게, 적은 진동을 원한다면 Ki 값을 적당히 낮게 하고 Kp Kd 값을 상대적으로 낮은 값() 으로 설정한다면 목적에 맞으면서도 안정적인 응답을 얻어 낼 수 있을 것이다.B. Frequency domain 으로의 해석각각 Ki = 20 Ki=70 일때의 Nyquist 선도와 Bode plot 을 구하면 다음과 같다.위 결과에 의하면 Ki=20 일때의 Pm < 0 이기 때문에 불안정함을 나타내게 된다 그렇기 때문에 Ki=20 일때 요구하던 small oscillation 을 유지하면서도 안정적이 되게 하기 위해 Ki값을 어느정도 올려보았다.Ki = 35 로 수정하였더니 Bode plot에서 Pm = 137 deg가 나와 Pm > 0을 만족시켜 안정된 응답이 도출됨을 알 수 있다.
2,3 DOF System ReportI. 2 DOF SYSTEM (Undamped free vibration)a. Analytical method-대수적 접근위와 같은 2 DOF system의 지배방정식은 다음과 같다.위의 지배방정식의 근을 밑의 식과 같다 가정하고m1,m2 와 k1,k2,k3 는 서로 선형화 해서 표현이 가능 하므로이렇게 나타낼 수 있다. 위의 두가지 수식을 지배방정식에 대입하여 풀어내면 다음과 같다.이를 vector [X]에 맞춰 모으고해서 풀어내면이때 [X]가 Nontrivial solution을 갖기 위해선 앞 행렬의 행렬식이 0이 되어야 한다. 이를 특성방정식 Characteristic equation 이라 한다. 하여튼 그렇게 λ에 대해 구하면 이를 통해 위의 식에서 제시된 형태로 [x] 를 표현 할 수 있고 Eigen vector 역시 구할 수 있다. 이때 위의 변수를 계속 문자로 놓고 풀면 식이 너무 복잡해져서 문제를 풀기 위한 상수를 제시하고 이에 관하여 풀겠다.-초기 조건 및 상수 제시위의 m,k 를 이용해 λ를 계산하자.λ가 다음과 같을때 우리는 각각 Eigen vector [X]1 [X]2 를 구해 낼 수 있다.위의 값을 이용하여 [x]를 표시하면상수 c1,c2,a1,a2 는 Initial condition 을 통해 구할 수 있다. 이를 통해 얻어진 식은 다음과 같다.-코딩Undamped 2DOF free vibration 의 m.file 중 Analytic method 의 부분이다. 손으로 계산한 식을 써서 graph로 그려내기까지의 과정이다.ti=0;tf=10;dt=0.1;t=ti:dt:tf;x1=-1.492*cos(2.05*t)*-0.860-(0.371)*(0.764)*cos(1.15*t);x2=-1.492*cos(2.05*t)*.510-(0.371)*(0.644)*cos(1.15*t);X(:,1)=x1;X(:,2)=x2;plot(t,x1,'s')hold onplot(t,x2,'h')gridxlabel('time')method-이론적 풀이 방법위의 지배방정식을 4th Runge Kutta method 를 사용하여 풀기 위해선 먼저 차수가 2인 미분방정식의 차수를 낮춰야 한다. 이를 위해서 함수 f를 이용해 x들을 각각 표현 하였다 그 풀이는 다음과 같다.이렇게 얻어진 f1,f2,f3,f4를 함수로 정의내리고 이를 ode 45에 적용하면 numerical solution 을 구할 수 있다.-코딩runge4 함수를 정의내리기 위해 만든 m.file 의 코딩 이다.function df=runge4(t,f);df= zeros(4,1)df(1)=f(3);df(2)=f(4);df(3)=-3*f(1)+2*f(2);df(4)=1*f(1)-2.5*f(2);endtwo dof 함수 중 ode45 를 이용해 풀어내고 graph를 그리는데 까지의 코딩이다.ti=0;tf=10;dt=0.1;t=ti:dt:tf;x0=[1,-1,0,0];[t y]=ode45('runge4',t,x0);plot(t,y(:,1),'v')hold onplot(t,y(:,2),'x')hold on-Graph▽ -> x1X -> x2C. 비교2가지 그래프를 동일 평면상에 그려보면 두 방법의 결과가 동일함을 확인 할 수 있다.실선 ->analytical solution▽,X->Numerical solutionII. 2 DOF SYSTEM (Damped free vibration)a. Analytical method-대수적 접근위와 같은 2 DOF system의 지배방정식은 다음과 같다위의 지배방정식의 근을 밑의 식과 같다 가정하고이를 집어넣어으로 바꿔 쓸 수 있다. 이때 각각의 상수들(m,c,k)는 각각 의 비로 표시되어 질 수 있다.이를 만족시키는 λ 값을 찾고 이를 이용하여 행렬의 Eigen vector 값을 구하면 우리는의 값을 시간에 관하여 표현할 수 있다. 문자를 그대로 쓰기엔 너무 복잡하여 모델링을 할 수치를 대입하여 풀이를 진행하도록 하겠다.-초기 조건 및 상수 제시위의 수치를 지배방정식 에 대입하여 λ 값을 구위해선 앞의 행렬식이 0이 되어야 한다. 이를 이용하여 λ값을 구하면이제 각각의 λ 를 위의 식 에 대입해서 4개의 Eigen vector 를 계산할 차례이다.이때 λ 1,2 λ 3,4 가 각각 켤레 복소수 이기 때문에 다음과 같은 꼴로 풀어 낼 수 있다.이때 Initial condition 을 이용하면 각 변수를 구할 수 있다. 변수를 구한 후 작성한 최종 식은 다음과 같다.-GraphX -> x1O -> x2B. Numerical method-이론적 풀이 방법Undamped 2DOF system에 사용한 방법과 유사한 방법을 사용한다. 위의 지배방정식을 함수 f 를 이용하여 차수를 낮춘 후에 각각 x들을 표현하겠다.Undamped 와 비슷 하긴 하지만 중간에 damper 로 인한 항이 추가 되었다. 위에서 구한 f들을 ode45 함수에 넣어줌으로서 쉽게 해를 구할 수 있다.-코딩2 DOF system 과 비슷한 방법을 사용한다. 먼저 runge4_1 함수를 정의내리기 위해 만든 m.file 의 코딩이다.function df=runge4_1(t,f);df= zeros(4,1)df(1)=f(3);df(2)=f(4);df(3)=-0.3*f(3)+0.2*f(4)-3*f(1)+2*f(2);df(4)=0.1*f(3)-0.25*f(4)+1*f(1)-2.5*f(2);endTwodofdamp m.file 중 ode45 를 이용해 풀어내고 graph를 그리기 위한 코딩이다.ti=0;tf=10;dt=0.1;t=ti:dt:tf;x0=[1,-1,0,0];[t y]=ode45('runge4_1',t,x0);plot(t,y(:,1),'b')hold onplot(t,y(:,2),'r')gridxlabel('time')ylabel('displacement')-GraphBlue->x1Red->x2C. 비교2가지 그래프를 동일 평면상에 그려보면 Damped system 역시 Undemped system 과 동일하게 두 방법의 결과가 동일함을 확인 할 수 있다.Line -> Numerical(Undamped free vibration)a. Analytical method-대수적 접근위와 같은 3 DOF system의 지배방정식은 다음과 같다.앞에서 한 방식(Undamped 2DOF System)과 같이 위의 지배방정식의 근을 밑의 식과 같다 가정하자.위 식을 지배방정식에 대입하고 연산을 해주면 다음과 같다.위 식에서 앞에서 사용한 방법대로 [X]가 Non trivial solution 이 되기 위해선 앞 행렬의 행렬식이 0 이 되어야 한다. 이를 이용해 ω 를 구할 수 있다. 이 이후부턴 문자로 풀이가 난해하기 때문에 상수를 제시하여 풀이 하겠다.-초기 조건 및 상수 제시위의 지배방정식에 상수를 대입하면 다음과 같다.λ가 다음과 같을때 우리는 각각 Eigen vector [X]1 [X]2 를 구해 낼 수 있다. 이때 λ 는 ω^2 인 값이다.위의 값을 이용하여 [x]를 표시하면상수 c1 c2 c3는 위 식의 전개와 Initial condition을 이용해 구할 수 있다. 상수를 모두 구한 식은 다음과 같다.-코딩Threedof.m file 중 analytic method 부분이다. Eigen vector 를 l1,l2,l3, 를 통해 나타내고, 손으로 푼 식을 적어서 변위를 표현하였다.ti=0;tf=10;dt=0.1;t=ti:dt:tf;l1=[-0.2114;-0.2507;-0.1720];l2=[0.3814;-0.0953;-0.1907];l3=[0.0993;-0.1675;0.3661];x=+0.59*l1*cos(0.7923*t)+1.907*l2*cos(1.871*t)+4.002*l3*cos(2.524*t)plot(t,x,'o')gridxlabel('time')ylabel('displacement')-GraphBlue->x1Red->x2Green->x3B. Numerical method-이론적 풀이 방법3 DOF system을 Numerical하게 푸는 방법은 위에 풀이한 2 DOF system과 크게 다르지 않다 다만 변위값이 2개에서 3개로 늘어. 위의 방법과 동일하게 f를 이용해 차수를 낮춰보는 방법은 다음과 같다.-코딩runge4_2함수를 정의내리기 위해 만든 m.file 의 코딩 이다.function df=runge4_2(t,f);df= zeros(6,1)df(1)=f(4);df(2)=f(5);df(3)=f(6);df(4)=-3*f(1)+2*f(2);df(5)=f(1)-2.5*f(2)+1.5*f(3);df(6)=3*f(2)-5*f(3);endThreedof m.file 중 ode45 를 이용해 풀어내고 graph를 그리기 위한 코딩이다.ti=0;tf=10;dt=0.1;t=ti:dt:tf;x0=[1,-1,1,0,0,0];[t y]=ode45('runge4_2',t,x0);plot(t,y(:,1),'b')hold onplot(t,y(:,2),'r')hold onplot(t,y(:,3),'g')hold ongridxlabel('time')ylabel('displacement')-GraphBlue -> x1Red -> x2Green -> x3C. 비교3 DOF System도 마찬가지로 두가지 결과가 일치함을 확인할 수 있다.Shape->analytical solutionLine ->Numerical solutionIV. Consideration이번 과제를 하면서 많이 어려웠던 점은 Analytic method를 푸는데 있어 굉장히 복잡한 계산이였다. 하면서 몇몇 새로운 기능(eig() 같은)을 익혀서 덜해지긴 했지만 상당히 많은 양의 계산이 필요하고, 사소한것에 실수를 자주해서 결과값이 안맞아 고생을 하기도 하였다. 이번 레포트를 통해 Runge Kutta method가 얼마나 강력한 풀이법인지 다시한번 알 수 있었다.3가지 종류의 진동에 대해 계산을 했는데 이중 가장 까다로운 종류는 damper가 달린 2DOF system 이였다. damper가 들어가는 순간 Eigen value 와 Eigen vector 구하기가 굉장히 힘들어 졌다. 차수가 더 높은 damping system 에 대하여서는 더 복.
Damped Viscous Free Vibration교수님:학번:학년:성명:제출일:1. 문제 풀이 과정a. Analytical method위의 미분방정식은 2차 선형 Homogeneous 미분방정식이다. 각 미지수와 그 미분한 값을 을 λ에 대한 식으로 바꾸어서 간단히 해결 할 수 있다. 이때 변환된 수식은 다음과 같다.이때 λ의 값이 a 라면으로 표현 할 수 있다. 각각 λ를 구하여 x에 대해 풀어 쓰면 다음과 같다.위의 수식 3을 이용하여 Analytical 한 해를 구할 수 있었다. λ가 중근 허근이 나올때 즉 ζ값이 1 이상일때 x값은 수식 4와 같다.(chapter 4 수업 print 참조)b. Numerical method (4th order Runge-Kutta method)수식1 의 미분방정식은 수식5의 2개의 1차 미분 방정식으로 변형되어 지고, 이는 4th order Runge Kutta method 를 이용해 쉽게 구할 수 있다. Runge Kutta method는 미분방정식을 푸는 하나의 방법으로써k값들에 가중치를 두어 평균을 구함으로서 구해낼 수 있다. Matlab 에서는 이를 직접 코딩할 수 있고 ODE45를 통해 풀어 낼 수 있는데 2가지 방법 모두 사용 하였다2. Sourcea. Hw1.m (본 함수 포함한 파일 코딩)function VTB2_6(m,k,dtype,dcoef,dt,tott,x0,v0);close allm=10;k=100;dtype=1;dcoef=1.5;dt=0.01;tott=10;x0=10;v0=0;syms y wn f wdn=floor(tott/dt+1);z(:,1)=[x0,v0]h=dt;t=0:dt:tott;k1=[0;0];k2=[0;0];k3=[0;0];k4=[0;0];wn=sqrt(k/m);wd=sqrt(abs(1-dcoef^2))*wn;if dcoef1y=exp(-dcoef*wn*t).*((x0*wn*(dcoef+sqrt(dcoef^2-1))+v0)/(2*wn*sqrt(dcoef^2-1))*cosh(sqrt(dcoef^2-1)*wn*t)+(-x0*wn*(dcoef-sqrt(dcoef^2-1))-v0)/(2*wn*sqrt(dcoef^2-1))*sinh(sqrt(dcoef^2-1)*wn*t))end%Analytical method 를 이용한 풀이ri=0.0;rf=10.0;rintval = [ri rf];bcinir = [x0,v0];[r x ] = ode45('vibrate', rintval, bcinir);%ODE 45 를 이용한 풀이if dtype==1for l1=1:(n-1);z1=z(:,l1);k1(1)=h*z1(2);k1(2)=h*(-wn^2*z1(1)-2*wn*dcoef*z1(2));k2(1)=h*(z1(2)+.5*k1(2));k2(2)=h*(-wn^2*(z1(1)+.5*k1(1))-2*wn*dcoef*(z1(2)+.5*k1(2)));k3(1)=h*(z1(2)+.5*k2(2));k3(2)=h*(-wn^2*(z1(1)+.5*k2(1))-2*wn*dcoef*(z1(2)+.5*k2(2)));k4(1)=h*(z1(2)+k3(2));k4(2)=h*(-wn^2*(z1(1)+k3(1))-2*wn*dcoef*(z1(2)+k3(2)));z(:,l1+1)=z(:,l1)+(k1+2*k2+2*k3+k4)/6;endendif dtype==2for l1=1:(n-1)z1=z(:,l1);k1(1)=h*(z1(2));k1(2)=h*(-wn^2*z1(1)-2*wn*dcoef*sign(z1(2)));k2(1)=h*(z1(2)+.5*k1(2));k2(2)=h*(-wn^2*(z1(1)+.5*k1(1)))-2*wn*dcoef*sign(z1(2)+.5*k1(2));k3(1)=h*(z1(2)+.5*k2(2));k3(2)=h*(-wn^2*(z1(2)+.5*(k2(1))-2*wn*dcoef*sign(z1(2)+.5*k2(2))));k4(1)=h*(z1(2)+k3(2));k4(2)=h*(-wn^2*(z1(1)+k3(1))-2*wn*dcoef*sign(z1(2)+k3(2)));z(:,l1+1)=z(:,l1)+1/6*(k1+2*k2+2*k3+k4);if abs(z1(2))k*z1(1)z(:,l1+1)=z(:,l1);z(2,l1+1)=0;endendendif dtype==3for l1=1:(n-1);z1=z(:,l1);k1(1)=h*(z1(2));k1(2)=h*(-wn^2*z1(1)-2*wn*dcoef*sign(z1(2))*(z1(2))^2);k2(1)=h*(z1(2)+.5*k1(2));k2(2)=h*(-wn^2*(z1(1)+.5*k1(1))-2*wn*dcoef*sign(z1(2)+.5*k1(2))*(z1(2)+.5*k1(2))^2);k3(1)=h*(z1(2)+.5*k2(2));k3(2)=h*(-wn^2*(z1(1)+.5*k2(1))-2*wn*dcoef*sign(z1(2)+.5*k2(2))*(z1(2)+.5*k2(2))^2);k4(1)=h*(z1(2)+k3(2));k4(2)=h*(-wn^2*(z1(1)+k3(1))-2*wn*dcoef*sign(z1(2)+k3(2))*(z1(2)+k3(2))^2);z(:,l1+1)=z(:,l1)+1/6*(k1+2*k2+2*k3+k4);endend%Runge Kutta method 를 직접 코딩한 풀이plot(t,y)figure(2)plot(t,z(1,:),'r')%figure(3)%plot(r,x(:,1),'g')%결과값 graph로 출렭b. Vibrate.m (ODE 45 를 사용하기 위해 함수를 정의하기 위한 m.file)이때 들어가 있는 상수 값은 ζ=0.3,ωn=10 인 상황이다. 이 함수속에 원래 변수로 하여 한번에 되도록 하려고 하였으나 매틀랩 실력이 부족하여 변수를 바꿀때마다 여기도 바꿔 줘야 했다.3. Conclusiona. 각 방법별 graphAnalytical methodNumerical method(coding)Numerical method(ode 45)m=10, k=1000, dtype=1, dcoef=0.3, dt=0.01, tott=10, x0=10, v0=0위의 표에서 보이는 그래프들과 같이 3가지 방법 모두 동일한 Graph 를 얻을 수 있다. 이는 dt가 충분히 작을 때 numerical 한 값이 analytical 한 값을 충분히 정확하게 근사 할 수 있음을 의미한다.b. Ω값의 변화에 따른 그래프.(m=10, k=100)Analytical methodNumerical methodΩ=0.3Ω=0.7Ω=1.0Ω=1.34. ConsiderationNumerical 한 방법은 Analytical 하게 접근하기 어려운 여러 비선형 문제를 쉽게 해결할 수 있게 해주는 훌륭한 수학적 도구이다. 이번 레포트에서는 선형 미분방정식의 Analytic 해와 Numerical 해를 비교해 봄으로써 그 사용 방법과 정확도를 알 수 있었다. 이때 주의해야 할 것 중 하나가 dt 값이었는데(시간간격) dt값이 충분히 조밀하지 못하면 무척 엉뚱한 값이 튀어나오기 때문에 주의해야한다. 밑의 4 개의 그래프로 dt값에 따른 오차를 확인 할 수 있다.dt=0.01 적절한 dt 값이 주어진 상태이다.dt=0.1 약간 큰 dt값이 주어진 상태이다. 그래프가 약간 날카로워졌다.dt=0.2 dt값이 임계치를 넘어간 경우이다. 진동에 대한 정확한 정보를 얻을 수 없다.dt=0.3 너무 큰 dt 값이 주어진 상태이다. 전혀 맞지 않는 모양의 그래프가 나온다.5.Reference-기계진동학 Balakumar Balachandran / Edward B. Magrab-응용수치 해석 Steven C. Chapra-Matlab application : Differential Equations and Boundary Value Problems-> http://people.rit.edu/pnveme/Matlab_Course/Matlab_App_ODE.html-Chap 4: S-DOF system Free Vibration (수업 ppt)
-모래언덕과 전쟁-‘세상은 생각보다 단순하다’ 를 읽고평평한 판 위에 모래알이 하나 둘 떨어진다. 어느새 높이 산을 이룬다. 산중턱에서 모래알이 살짝 휩쓸린다. 한번 두 번 그렇게 크고 작은 휩쓸림이 계속 이뤄지다가 단 하나의 모래알이 더 떨어졌을 때 엄청나게 큰 산사태가 벌어진다. 우리는 이때 마지막에 떨어진 단 하나의 모래알이 이렇게 거대한 산사태를 만들어 낼지 알 수 있었을까? 만약 이에 대해 예라고 대답할 수 있다면 우리는 인류 최대의 비극인 1차 세계대전과 2차 세계대전을 막을 수 있었을 것이며 87년도에 일어난 증시사상 최악의 폭락사태나 2009년 ‘서브프라임모기지‘로 인한 세계 증시의 붕괴 역시 예측할 수 있었을 것이다. 왜 우리는 이러한 역사의 극적인 순간을 이전의 순환되는 정보로부터 예측할 수 없는 것일까? 이 책의 필자는 이 물음에 대해 현상을 극단적으로 단순화한 모델을 이용한 게임을 통해 우리에게 이러한 현상들 속에서 하나의 유사성을 보여주고 있다.모래알이 떨어지는 각각의 사건은 이전이나 지금이나 다를 바 없다. 하지만 중요한 사실은 언덕의 상태, 즉 언덕이 임계상태인지 아닌지, 모래알이 떨어지기 직전까지 언덕이 받아온 스트레스는 어느 정도 인지이다. 임계상태의 언덕은 각각 모래알이 떨어질 때 그 위치에 따라 단지 작은 휩쓸림만으로 끝날 수도 있고 백만 개가 넘는 알갱이가 무너져 버릴 수도 있다 . 이러한 결과는 무너져 내리는 알갱이의 크기와 발생 빈도의 관계에서도 확인 할 수 있다. 무너져 내리는 알갱이가 2배로 많아지면 발생 빈도는 2배쯤 감소하게 된다. 이러한 관계를 멱함수 관계라고 한다. 이러한 관계를 갖는 현상의 두드러진 특징은 현상의 전형적인 크기가 존재하지 않는다는 것이다. 이는 다음번 떨어질 모래알이 겨우 몇 개의 모래알만 떨어뜨릴 수도 있고 혹은 엄청나게 큰 산사태의 시발점이 될 수도 있음을 의미한다.모래알과 언덕 붕괴사이의 관계는 자연계의 많은 부분을 설명할 수 있는 열쇠가 된다. 가장 유명한 사례 중 하나는 지진의 세기와 발생 빈도에 따른 구텐베르크 - 리히터의 법칙 이다. 이는 지진의 세기가 2배 강해질수록 발생빈도는 4배 줄어든다는 법칙이다. 지금까지의 지진예측을 전부 빗나가게 한 원인도 바로 이 성질 때문이다. 모래알이 끊임없이 떨어지듯 맨틀은 끊임없이 대류 하는데 이로 인해 지각이 단 수 밀리미터만 움직일 수도 아니면 95년의 고베지진과 같은 대지진이 일어날 수도 있다. 원인은 자극의 세기나 이전 지진의 발생경향 따위가 아니다. 여태까지 조그마한 지진과 지각의 움직임으로 인해 쌓인 스트레스, 바로 그 스트레스 지점의 자극이 커다란 산사태, 대지진의 가장 직접적인 원인이라고 할 수 있다. 혹자는 그럼 그러한 지점에 자극이 가해질 확률을 구하면 되지 않겠냐고 반문할 수 있을 것이다. 하지만 이러한 스트레스 지점이 아닌 다른 곳에 가해진 자극은 결국 다른 지점의 스트레스를 쌓을 것이고 이러한 단순 회피는 연쇄작용으로 인한 더욱 큰 재앙만을 부를 뿐이다. 그 뿐만 아니라 이러한 지점을 결정하는 요인은 굉장히 다양한 요인이 시간에 따라 복잡하게 엉켜있기 때문에 이러한 위치를 정확하게 예측하는 것은 불가능에 가깝다.지진뿐만 아니라 자연계 거의 대부분의 복잡 계는 모래알과 언덕의 관계를 가지고 있다. 자석 내 원자스핀의 쏠림 정도와 그 빈도, 생물 진화속도와 그 빈도, 개체 수 급증과 그 빈도, 종의 멸종과 빈도 역시 모두 멱함수 관계를 가지고 있다. 이들의 공통점은 개체들이 단순히 홀로 무언가에 작용하는 것이 아닌 서로간의 상호작용을 통해 연쇄반응을 일으킨다는 것이다. 하나의 종의 진화, 급증, 멸종은 그물처럼 총총히 짜인 자연계 대부분에 영향을 미치며 이는 마치 임계상태의 모래언덕처럼 꿈쩍도 안할 수도, 산사태처럼 모조리 박살날 수 도 있다.과연 인류의 역사에선 이들과의 유사성을 찾을 수 없을까? 인간은 자유의지를 가지고 있기 때문에 무작위로 가해지는 자극에 능동적으로 대체하여 최악의 결과를 막을 수 있는 것처럼 보이지만 사람들을 모아놓은 집단에서는 그 얘기가 달라진다. 주식시장의 폭락, 상승 빈도와 혁명의 발생빈도, 전쟁의 규모 빈도의 관계까지 모래알과 언덕의 관계와 유사하다. 인류사회의 역사들은 얼핏 보면 굉장히 중요한 원인에 의해서 위대하고 특별한 사람들에 의해 쓰여 진 것처럼 보이지만 실제론 적절한 시기와 적절한 위치에 우발적으로 일어난 사건들에 의해 혁명이 일어나고 세계대전이 일어나고 주가가 폭락하고 역사의 큰 물줄기가 바뀌게 된다. 마치 임계상태의 모래언덕처럼, 자그마한 모래 한 알에도 모래언덕 전체의 모양이 바뀌게 돼버린다. 이를 이해하기 위해선 모래알이 떨어짐에만 집중하면 안 되는 것처럼 역사의 물줄기를 바꿔놓은 사건들을 정확히 이해하기 위해선 사건의 원인에만 집중하는 것이 아닌 규칙성 없어 보이는 사건들의 연결성과 그 유사성에 대해 탐구해봐야 할 것이다.