diff --git a/AERO3220_HW4.m b/AERO3220_HW4.m new file mode 100644 index 0000000..d396f2a --- /dev/null +++ b/AERO3220_HW4.m @@ -0,0 +1,207 @@ +% Hayden Molloy +% AERO3220 HW4 +% Due 18 Feb 2026 + +clear; close all; clc; + + +%% Initial Parameters + + k1 = 100; % N/m + k2 = 50; % N/m + m1 = 500; % kg + m2 = 250; % kg + +% Time Parameters + t0 = 0; + tf = 200; + dt = 0.01; + tspan = t0:dt:tf; + +% Initial conditions [z1 z1dot z2 z2dot] + x0 = [0; 0; 0; 0]; + +% Initialize Force Functions + F_1 = @(t) 100; % N + F_2 = @(t) 100*sin(0.25*t); % N + F_3 = @(t) 100*sin(0.50*t); % N + + +%% Define Cases and ODEs + +% CASE 1: + b = 100; + F = F_1; +[t1,x1] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0); + +% Case 2: + b = 100; + F = F_2; +[t2,x2] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0); + +% CASE 3: + b = 0; + F = F_2; +[t3,x3] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0); + +% Case 4: + b = 100; + F = F_3; +[t4,x4] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0); + +% Case 5: + b = 0; + F = F_3; +[t5,x5] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0); + + +%% Report Results + +% Create Legend +lgnd = { ... + 'Case 1: F=100, b=100', ... + 'Case 2: F=100sin(0.25t), b=100', ... + 'Case 3: F=100sin(0.25t), b=0', ... + 'Case 4: F=100sin(0.50t), b=100', ... + 'Case 5: F=100sin(0.50t), b=0'}; + + +% Plot 1: Top mass displacement vs. Time +figure(1); + +plot(t1,x1(:,1),'LineWidth',1.3); +hold on; +plot(t2,x2(:,1),'LineWidth',1.3); +plot(t3,x3(:,1),'LineWidth',1.3); +plot(t4,x4(:,1),'LineWidth',1.3); +plot(t5,x5(:,1),'LineWidth',1.3); +grid on; +xlabel('Time (s)'); +ylabel('z_1 (m)'); +title('Top Mass Displacement z_1(t)'); +legend(lgnd,'Location','southoutside'); + + +% Plot 2: Bottom mass displacement vs. Time +figure(2); + +plot(t1,x1(:,3),'LineWidth',1.3); +hold on; +plot(t2,x2(:,3),'LineWidth',1.3); + +plot(t3,x3(:,3),'LineWidth',1.3); + +plot(t4,x4(:,3),'LineWidth',1.3); + +plot(t5,x5(:,3),'LineWidth',1.3); + +grid on; +xlabel('Time (s)'); +ylabel('z_2 (m)'); +title('Bottom Mass Displacement z_2(t)'); +legend(lgnd,'Location','southoutside'); + + +% Plot 3: Top mass velocity vs. Time +figure(3); + +plot(t1,x1(:,2),'LineWidth',1.3); +hold on; +plot(t2,x2(:,2),'LineWidth',1.3); + +plot(t3,x3(:,2),'LineWidth',1.3); + +plot(t4,x4(:,2),'LineWidth',1.3); + +plot(t5,x5(:,2),'LineWidth',1.3); + +grid on; +xlabel('Time (s)'); +ylabel('dz_1/dt (m/s)'); +title('Top Mass Velocity z_1'''); +legend(lgnd,'Location','southoutside'); + + +% PLOT 4: Bottom mass velocity vs. Time +figure(4); + +plot(t1,x1(:,4),'LineWidth',1.3); +hold on; +plot(t2,x2(:,4),'LineWidth',1.3); + +plot(t3,x3(:,4),'LineWidth',1.3); + +plot(t4,x4(:,4),'LineWidth',1.3); + +plot(t5,x5(:,4),'LineWidth',1.3); + +grid on; +xlabel('Time (s)'); +ylabel('dz_2/dt (m/s)'); +title('Bottom Mass Velocity z_2'''); +legend(lgnd,'Location','southoutside'); + + +% PLOT 5: z1(t) and z2(t) +% (solid = z1, dashed = z2) + +c1 = [0 114 178]/255; +c2 = [213 94 0]/255; +c3 = [230 159 0]/255; +c4 = [117 112 179]/255; +c5 = [0 158 115]/255; + +figure(5); + +plot(t1,x1(:,1),'-','LineWidth',1.2, 'Color', c1); +hold on; +plot(t1,x1(:,3),'--','LineWidth',1.2, 'Color', c1); + +plot(t2,x2(:,1),'-','LineWidth',1.2, 'Color', c2); +plot(t2,x2(:,3),'--','LineWidth',1.2, 'Color', c2); + +plot(t3,x3(:,1),'-','LineWidth',1.2, 'Color', c3); +plot(t3,x3(:,3),'--','LineWidth',1.2, 'Color', c3); + +plot(t4,x4(:,1),'-','LineWidth',1.2, 'Color', c4); +plot(t4,x4(:,3),'--','LineWidth',1.2, 'Color', c4); + +plot(t5,x5(:,1),'-','LineWidth',1.2, 'Color', c5); +plot(t5,x5(:,3),'--','LineWidth',1.2, 'Color', c5); + +grid on; +xlabel('Time (s)'); +ylabel('Displacement (m)'); +title('Top and Bottom Mass Displacements'); + +legend({ ... + 'Case 1 z_1', ... + 'Case 1 z_2', ... + 'Case 2 z_1', ... + 'Case 2 z_2', ... + 'Case 3 z_1', ... + 'Case 3 z_2', ... + 'Case 4 z_1', ... + 'Case 4 z_2', ... + 'Case 5 z_1', ... + 'Case 5 z_2'}, ... + 'Location','eastoutside'); + + +%% Local ODE function +function dx = eom_two_mass(t, x, m1, m2, k1, k2, b, Ffun) + z1 = x(1); + z1dot = x(2); + z2 = x(3); + z2dot = x(4); + + u = Ffun(t); + + % m1*z1ddot + k1*z1 - k1*z2 = 0 + z1ddot = (-k1/m1)*z1 + (k1/m1)*z2; + + % m2*z2ddot + b*z2dot + (k1+k2)*z2 - k1*z1 = u + z2ddot = (k1/m2)*z1 - ((k1+k2)/m2)*z2 - (b/m2)*z2dot + (1/m2)*u; + + dx = [z1dot; z1ddot; z2dot; z2ddot]; +end \ No newline at end of file diff --git a/AERO3220_HW4.pdf b/AERO3220_HW4.pdf new file mode 100644 index 0000000..53eaee7 Binary files /dev/null and b/AERO3220_HW4.pdf differ diff --git a/AERO3220_HW4_block_diagram.md b/AERO3220_HW4_block_diagram.md new file mode 100644 index 0000000..bdaa343 --- /dev/null +++ b/AERO3220_HW4_block_diagram.md @@ -0,0 +1,90 @@ +# AERO3220 HW4 Program Flow Block Diagram + +```mermaid +flowchart TD + A([Start Script]) --> B[Clear Workspace: clear, close all, clc] + B --> C[Set Parameters: k1, k2, m1, m2] + C --> D[Set Time Grid: t0, tf, dt, tspan] + D --> E[Set Initial State: x0 = [z1 z1dot z2 z2dot]^T] + E --> F[Define Input Force Functions: F1, F2, F3] + + F --> G[Case 1: b=100, F=F1] + G --> H[Run ode45 with eom_two_mass] + + H --> I[Case 2: b=100, F=F2] + I --> J[Run ode45 with eom_two_mass] + + J --> K[Case 3: b=0, F=F2] + K --> L[Run ode45 with eom_two_mass] + + L --> M[Case 4: b=100, F=F3] + M --> N[Run ode45 with eom_two_mass] + + N --> O[Case 5: b=0, F=F3] + O --> P[Run ode45 with eom_two_mass] + + P --> Q[Build Legend Strings] + Q --> R[Plot Figure 1: z1 vs time (5 cases)] + R --> S[Plot Figure 2: z2 vs time (5 cases)] + S --> T[Plot Figure 3: z1dot vs time (5 cases)] + T --> U[Plot Figure 4: z2dot vs time (5 cases)] + U --> V[Plot Figure 5: z1 and z2 together by case] + + V --> W[Local Function: eom_two_mass(t,x,...)] + W --> X[Unpack States: z1, z1dot, z2, z2dot] + X --> Y[Evaluate Input: u = Ffun(t)] + Y --> Z[Compute Accelerations: z1ddot and z2ddot] + Z --> AA[Return State Derivative: dx = [z1dot z1ddot z2dot z2ddot]^T] + AA --> AB([End]) +``` + +## ODE Model Used in Each Case + +\[ +\ddot{z}_1 = -\frac{k_1}{m_1}z_1 + \frac{k_1}{m_1}z_2 +\] + +\[ +\ddot{z}_2 = \frac{k_1}{m_2}z_1 - \frac{k_1+k_2}{m_2}z_2 - \frac{b}{m_2}\dot{z}_2 + \frac{1}{m_2}u(t) +\] + +where \(u(t)\) is one of the three forcing functions depending on the case. + +## Assignment-Style Block Diagram (Recommended For Submission) + +```mermaid +flowchart TD + A([Start]) --> B[Define constants and masses: k1, k2, m1, m2] + B --> C[Define simulation time: t0, tf, dt, tspan] + C --> D[Set initial state vector: x0 = [z1 z1dot z2 z2dot]^T] + D --> E[Define force inputs: F1=100, F2=100sin(0.25t), F3=100sin(0.50t)] + + E --> F[Create case table: (Case, b, Force)] + F --> G{i <= Number of Cases?} + + G -- Yes --> H[Load case i values: b_i, F_i] + H --> I[Call ode45 with eom_two_mass and current case] + I --> J[Store outputs: t_i and x_i] + J --> K[i = i + 1] + K --> G + + G -- No --> L[Generate legend labels] + L --> M[Plot z1 vs time for all cases] + M --> N[Plot z2 vs time for all cases] + N --> O[Plot z1dot vs time for all cases] + O --> P[Plot z2dot vs time for all cases] + P --> Q[Plot combined z1 and z2 by case] + Q --> R([End]) +``` + +### Solver Subsystem (Inside `eom_two_mass`) + +```mermaid +flowchart LR + A1[Inputs: t, x, m1, m2, k1, k2, b, Ffun] --> A2[Unpack states: z1, z1dot, z2, z2dot] + A2 --> A3[Compute input force: u = Ffun(t)] + A3 --> A4[Compute z1ddot from mass 1 equation] + A4 --> A5[Compute z2ddot from mass 2 equation] + A5 --> A6[Assemble derivative: dx = [z1dot z1ddot z2dot z2ddot]^T] + A6 --> A7[Return dx to ode45] +```