Question

An engineering company is designing and testing a car suspension system. The system has a conventional suspension design, consisting of a shock-absorber and spring at each wheel. The shock-absorber provides a damping effect that is proportional to the vertical speed (i.e. up and down motion) of the wheel, and the spring's resistance is proportional to the vertical displacement of the wheel.

The design team analyses the suspension system in two separate parts: Part 1 - dynamics of the spring-damper system; Part 2 - stress analysis of the coil spring.

PART 1

The company would like to investigate how the actual performance of the suspension system compares to that of theory. Dynamic theory states that the equation of motion of the wheel in the vertical direction can be expressed as,

  $LaTeX: m\frac{d^2y}{dt^2}=-c\frac{dy}{dt}-ky$ or $LaTeX: \frac{d^2y}{dt^2}+\frac{c}{m}\frac{dy}{dt}+\frac{k}{m}y=0$

where m is the mass of the wheel, k is the spring stiffness and c is the damping coefficient.

This is a second order differential equation, and if for example, the car hits a hole at t = 0, such that it is displaced from its equilibrium position with y = y₀, and dy/dt = 0, it will have a solution of the form,

where    and $LaTeX: n=\frac{c}{2m}$ provided   $LaTeX: \frac{k}{m}>\frac{c^2}{4m^2}$

To analyse the actual performance of the design, the suspension system was built and tested with the following displacements of the wheel recorded during the time period of 2.6 and 3.1 seconds (this dataset is known to contain experimental error):

Time (s)	2.6	2.65	2.7	2.75	2.8	2.85	2.9	2.95	3	3.05	3.1
Displacement (mm)	0.074	0.094	0.106	0.113	0.127	0.145	0.145	0.148	0.154	0.162	0.158

The design team are interested if the dataset can be characterised by expressions which are less complex that the above theory.

PART 2

The design team performed stress analysis on the coil spring, and it was determined that the stress at a specific point could be characterised by:

Similar calculations were performed to calculate the maximum principal stress at 250 locations along a non-linear path through the coil spring. The following frequency distribution of the calculated maximum principal stresses was produced:

Max. Principal Stress (MPa)	0.5	1	1.5	2	2.5	3	3.5	4	4.5	5
Frequency	10	15	19	24	36	50	42	29	20	5

The design team is interested in the significance of these results and therefore need to calculate the area under the graph which contains this dataset.

The design team believe any further calculations of principal stresses and determining the area under the graph will be time consuming (and potentially error prone), thus would like an automated method for performing this.

Submission

Create two separate MATLAB scripts for Part 1 and Part 2:

PART 1 - Create a MATLAB script which is capable of performing the following:

• Plotting the theoretical displacement of the wheel which shows the first 6 roots when the suspension system has the following specifications:
  • mass acting on each wheel = 3.6 x 10³ kg
  • c = 1.5 x 10³ Ns/m
  • k = 1.5 x 10⁴ N/m                          
  • initial displacement = 0.3 m
• Calculating the time for the first 3 occasions the wheel passes through the equilibrium position (i.e. the root);
• Plotting a graph of the experimental dataset of the wheel displacement;
• Calculating the value of the coefficient of multiple determination ( R 2 ) and the standard error ($LaTeX: \sigma_E$ σ E ) when presuming the experimental dataset is:
  • Linear (i.e. y = a + b x )
  • Characterised by a saturation growth equation (i.e. $LaTeX: y=\frac{ax}{b+x}$ y = a x b + x )

PART 2 – Create a MATLAB script that is capable for performing the following:

• Calculating the maximum principal stress at the specific point shown;
• Plot the frequency distribution of the calculated maximum principal stresses and fit a natural spline through the dataset;
• Calculate the area under the natural spline (between the data points).

Make sure all sections of the MATLAB script are annotated to explain the role of each line within the script.

v(t) -(o cosin(t) 4m2

Answer 1

PART --- 1

Code :=

%%%%%%%%%% Theroritical Ploting and Roots %%%%%%%%%%%%%
m = 3.6*10^3; c = 1.5*10^3;   k = 1.5*10^4; y0 = 0.3;
p = sqrt((k/m)-(c^2/(4*m^2)));
n = c/(2*m);
y = @(t) exp(-n*t).*(y0*cos(p*t) + y0*(n/p)*sin(p*t));

t = linspace(0,15);
plot(t,y(t))    % ploting theritical displacement
yline(0);       % Ploting zero line to identify root locations
xlabel("Time [s]"); ylabel("Displacement [m]"); title("Theoritical Plot")
% Form figure it's evident that first 3 roots are near 0.5, 2 and 4
% To calculate exact values we will use fzero function
root1 = fzero(y,0.5);   % 1st root using initial guess of 0.5
root2 = fzero(y,2);     % 2nd root using initial guess of 2
root3 = fzero(y,4);     % 3rd root using initial guess of 4
fprintf("\n First 3 roots : [%.4f %.4f %.4f] sec",[root1,root2,root3])
%%%%%%%%% Experimental Calculations %%%%%%%%%%%%%%%
warning off
Time = [2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3 3.05 3.1];
Disp = [0.074 0.094 0.106 0.113 0.127 0.145 0.145 0.148 0.154 0.162 0.158];
figure(2)
plot(Time,Disp,'*-')

[lin_Fit,gof_LinFit] = fit(Time',Disp','poly1');    % Linear Fit

% defininig fit type for saturation growth
myfittype = fittype('a*x/(b+x)','coefficients',{'a','b'});
[sat_growth_fit,gof_SatFit] = fit(Time',Disp',myfittype);
hold on
time_plot = linspace(2.6,3.1);
plot(time_plot,lin_Fit(time_plot),'LineWidth',1.2);     % Ploting Linear Fit
plot(time_plot,sat_growth_fit(time_plot),'LineWidth',1.2);   % ploting saturation growth fit
xlabel("Time"); ylabel("Displacement");
legend(["Experimental Data","Linear Fit","Saturation Growth fit"],'Location','NorthWest')
title("Experimental Plot");

% R2 and standard error for both
R2_LinFit = gof_LinFit.rsquare;     % Extracting Rsquare from goodness of fit (gof)
Std_err_LinFit = gof_LinFit.rmse;   % Extracting Standard error from gof

R2_SatFit = gof_SatFit.rsquare;
Std_err_SatFit = gof_SatFit.rmse;

fprintf("\n\n For Linear Fit : R2 = %.4f , Standard Error = %.4f", R2_LinFit,Std_err_LinFit);
fprintf("\n For Saturation Growth Fit : R2 = %.4f , Standard Error = %.4f\n", R2_SatFit,Std_err_SatFit);
warning on

Results:-

First 3 roots : [0.8239 2.3711 3.91821 sec F r Linear Fit : R2 = 0.9172 , Standard Error = o.ooee For Saturation Growth Fit : R2 = 0.4295 , Standard Error = 0.0232

Theoritical Plot 0.3 0.2 E 0.1 0 -0.1 -0.2 0.3 0 5 10 15 Time [s]

Experimental Plot 0.18 Experimental Data Linear Fit Saturation Growth fit 0.16 0.14 0.12 0.1 0.08 0.06 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3 3.05 3.1 Time

PART----- 2

Max Principal Stress for given condition : 5.3652 MPa

sigma_y3.1; sigma-z = 4.5; tau_yz - -1.4; % Formula for maximum principal Stress fprintf("\n Max Principal Stre 33 for given conditions %.4f MPa", sigma-max-P) : igma 0.5:0.5:5; req- [10 15 19 24 36 50 42 29 20 5]: plot (sigma, freg) % Ploting Frequncy Distribution pline_fit - fit (sigma',freq',"SmoothingSpline'): hold on p-plot (spline_fit,sigma, freq) p(1) .MarkerSize 15; p (2) . Linewidth = 1.6; legend (t"Distribution","Data Points", "Fitted Spline", 'Location','Best' xlabel ("sigma [MPa]") ylabel ("Frequency":grid orn % Fitting spline t data % Ploting Fitted Spline For area under curve we will make more fine data to get accurate value sigma_finer - linspace (0.5, 5,100) freq-finer-spline-fit (sigma-finer); % Finding Corresponding frequencies by fitted spline % generating 100 points for 31gma between 0.5 & 5 trap2(sigma-finer, freq-finer); % Finding area under curve using trapz area - fprintf("\n Area under Curve %.4f", area) % Priniting area -

Area under Curve = 121.9053

50 Distribution Data Points 45 _Fitted Spline 40 35 C 30 25 20 15 10 0.5 1 1.5 22.5 3 3.5 44.55 5.5 ơ [MPal

Text :-

sigma_y = 3.1;
sigma_z = 4.5;
tau_yz = -1.4;
% Formula for maximum principal Stress
sigma_max_P = (sigma_y + sigma_z)/2 + sqrt(((sigma_y - sigma_z)/2)^2 + tau_yz^2);

sigma = 0.5:0.5:5;
freq = [10 15 19 24 36 50 42 29 20 5];

plot(sigma,freq)        % Ploting Frequncy Distribution

spline_fit = fit(sigma',freq','SmoothingSpline');       % Fitting spline to data
hold on
p = plot(spline_fit,sigma,freq);        % Ploting Fitted Spline
p(1).MarkerSize = 15;
p(2).LineWidth = 1.6;
legend(["Distribution","Data Points","Fitted Spline"],'Location','Best')
xlabel("\sigma [MPa]");     ylabel("Frequency");    grid on

% For area under curve we will make more fine data to get accurate value
sigma_finer = linspace(0.5,5,100);      % generating 100 points for sigma between 0.5 & 5
freq_finer = spline_fit(sigma_finer);   % Finding Corresponding frequencies by fitted spline

area = trapz(sigma_finer,freq_finer);   % Finding area under curve using trapz

fprintf("\n Area under Curve = %.4f",area)  % Priniting area