이 문서는 MATLAB 코드를 올려놓은 자료입니다.이 코드로 해당 강의에서 전체 1등을 하였고, 가장 기본적인 형태의 스털링 엔진의 movement를 구현하였습니다.아마 개인 term project 정도로 내시기에 적당한 분량이라 생각 되고 혹시 team project라면 여기에 4절링크 구조 하나 정도 더 추가하시는 것을 추천 드립니다. 기본적인 형태는 식만 추가하면 비슷할 것 같습니다.해당하는 시간, 초기 (각)속도 및 (각)가속도를 입력 시에 4개의 그래프로 movement, 변위 또는 각/(각)속도/(각)가속도를 표현하였습니다. 또한 동영상 형태로 해당 그래프의 모습이 저장되기도 합니다.미리보기 방지를 위해 이렇게 글을 적게 되었습니다. 감사합니다.clc; clear all;% 링크 길이a=0.065; b=0.19;c=0.16; r=0.0475; % r : 실린더, 피스톤 반지름h=0.21; l=0.3; % h : 실린더 높이% l : 밑면에서 크랭크 중심 사이 거리t=0; step=0.05;Vol_min=100; Vol_max=0; % Vol : 행정 부피t_min=0; t_max=0; % 행정 부피가 최소, 최대일 때의 시간t_=input('작동 시간을 입력해주세요.[s]: ');omega2_initial=input('각속도를 입력해주세요.[rad/s]: ');omega2=omega2_initial;alpha2=input('각가속도를 입력해주세요.[rad/s^2]: ');TH2_initial=pi/2; % 90도에서 운동시작i=0;fprintf('n 1: 위치, 선속도, 선가속도 2: 각도, 각속도, 각가속도n')CASE=input('무엇을 보시겠습니까?: ');while(CASE ~=1 && CASE ~=2)fprintf('잘못입력하였습니다. 다시 입력해주세요.n')CASE=input('무엇을 보시겠습니까?: ');end% data set initialtime=zeros();data_d=zeros();data_f=zeros();data_TH3=zeros();data_TH4=zeros();data_omega3=zeros();data_omega4=zeros();data_alpha3=zeros();data_alpha4=zeros();data_d_velocity=zeros();data_f_velocity=zeros();data_d_accel=zeros();data_f_accel=zeros();det=zeros();v=VideoWriter('Stirling_Engine_MATLAB.avi'); % 비디오 입력open(v);clf; % 기존 figure 창 지우기while(1)i=i+1;if t>t_break;end% 실제 작동현황을 보여줌subplot(2,2,1);title('Stirling Engine')axis([-0.2 0.4 -0.2 0.4])xlabel('x(m)');ylabel('y(m)');grid on;cla;% loop inputTH2=TH2_initial+omega2_initial*t+0.5*alpha2*t^2;omega2=omega2_initial+alpha2*t;%% Theta 3,4TH3=2*pi-acos(a*cos(TH2)/b); % 이론 계산으로 구한 값은 예각인데, 실제로 구한 각은% 그림에서 보듯 180도가 넘는다. 즉, 2pi-theta 꼴이다.TH4=pi-asin(a*sin(TH2)/c); % 이론 걔산으로 구한 값은 예각인데, 실제로 구한 각은% 그림에서 보듯 둔각이다. 즉, pi-theta 꼴이다.% 각각의 각에 대해 sin, cos에 대입시 부호는 달라질 수 있지만, 크기는 같다.%% Distanced=a*sin(TH2)-b*sin(TH3);f=a*cos(TH2)-c*cos(TH4);%% Angular Velocityomega3=(a/b)*(sin(TH2)/sin(TH3))*omega2;omega4=(a/c)*(cos(TH2)/cos(TH4))*omega2;%% Velocityd_vel=a*omega2*cos(TH2)-b*omega3*cos(TH3);f_vel=-a*omega2*sin(TH2)+c*omega4*sin(TH4);%% Angular Accelerationalpha3=(a*alpha2*sin(TH2)+a*omega2^2*cos(TH2)-b*omega3^2*sin(TH3))/(b*sin(TH3));alpha4=(a*alpha2*cos(TH2)-a*omega2^2*sin(TH2)+c*omega4^2*sin(TH4))/(c*cos(TH4));%% Accelerationd_accel=a*alpha2*cos(TH2)-a*omega2^2*sin(TH2)-b*alpha3*cos(TH3)+b*omega3^2*sin(TH3);f_accel=-a*alpha2*sin(TH2)-a*omega2^2*cos(TH2)+c*alpha4*sin(TH4)+c*omega4^2*cos(TH4);%% Coordinate of the pointsx1=0; y1=0;x2=x1+a*cos(TH2); y2=y1+a*sin(TH2);x3=x1; y3=y1+d;x4=x1+f; y4=y1;x5=x1; y5=y1+l;x6=x1+l; y6=y1;%% Link by two pointsL1x=[x1 x2]; L1y=[y1 y2];L2x=[x2 x3]; L2y=[y2 y3];L3x=[x2 x4]; L3y=[y2 y4];%% Crank, CylinderL4x=[x3-r x3+r]; L4y=[y3 y3]; % 피스톤1L5x=[x4 x4]; L5y=[y4-r y4+r]; % 피스톤2L6x=[x5-r x5+r]; L6y=[y5 y5]; % 실린더1 밑면L7x=[x6 x6]; L7y=[y6+r y6-r]; % 실린더2 밑면L8x=[x5-r x5-r]; L8y=[y5-h y5]; % 실린더1 외벽L9x=[x5+r x5+r]; L9y=[y5-h y5]; % 실린더1 외벽L10x=[x6-h x6]; L10y=[y6-r y6-r]; % 실린더2 외벽L11x=[x6-h x6]; L11y=[y6+r y6+r]; % 실린더2 외벽%% Drawing of links, crank, cylinderline(L1x,L1y,'color','green','LineWidth',2)line(L2x,L2y,'color','red','LineWidth',2)line(L3x,L3y,'color','blue','LineWidth',2)line(L4x,L4y,'color','red','LineWidth',3)line(L5x,L5y,'color','blue','LineWidth',3)line(L6x,L6y,'color','black')line(L7x,L7y,'color','black')line(L8x,L8y,'color','black')line(L9x,L9y,'color','black')line(L10x,L10y,'color','black')line(L11x,L11y,'color','black')%% Data savedata_d(i)=y3;data_f(i)=x4;data_TH3(i)=TH3;data_TH4(i)=TH4;data_omega3(i)=omega3;data_omega4(i)=omega4;data_alpha3(i)=alpha3;data_alpha4(i)=alpha4;data_d_velocity(i)=d_vel;data_f_velocity(i)=f_vel;data_d_accel(i)=d_accel;data_f_accel(i)=f_accel;time(i)=t;t=t+step;switch CASEcase 1% 위치 해석subplot(2,2,2);plotresult(t_,time,data_d,data_f,...'time[s]','distance[m]')% 속도 해석subplot(2,2,3);plotresult(t_,time,data_d_velocity,data_f_velocity,...'time[s]','velocity[m/s]')% 가속도 해석subplot(2,2,4);plotresult(t_,time,data_d_accel,data_f_accel,...'time[s]','acceleration[m/s^2]')case 2% 각 해석subplot(2,2,2);plotresult(t_,time,data_TH3,data_TH4,...'time[s]','theta[rad]')% 각속도 해석subplot(2,2,3);plotresult(t_,time,data_omega3,data_omega4,...'time[s]','angular velocity[rad/s]')% 각가속도 해석subplot(2,2,4);plotresult(t_,time,data_alpha3,data_alpha4,...'time[s]','angular acceleration[rad/s^2]')otherwisebreak;endV=pi*(r^2)*((0.3-d)+(0.3-f));if Vol_min>VVol_min=V;t_min=t;endif Vol_max