\documentclass[journal, letterpaper]{IEEEtran}
\usepackage{graphicx}
\usepackage{url}
\usepackage{amsmath}
\usepackage{textgreek} % Greek to me, dawg
\usepackage{listings}
\usepackage{csvsimple}
\usepackage{longtable}
\begin{document}
\title{Automobile Suspension Design}
\author{Ian Arlen, Paul Gu\'enette, \& Ariel Levy}
\markboth{MAE 315 Dynamics of Machines}{}
\maketitle
\begin{abstract}
Vibration control is crucially important in ensuring a smooth ride for vehicle passengers. This study sought to design a suspension system for a car such that its mode of vibration would be predominantly bouncing at lower speeds, and primarily pitching at higher speeds. Our study used analytical and numerical methods to choose appropriate springs and dampers for the front and rear suspension. After an initial miscalculation, we succeeded in arriving at appropriate shocks for the vehicle with the desired modes of vibration at the specified frequencies.
We then assessed the maximum bouncing and pitching that the vehicle would experience under a specific set of conditions: travel at 40 km/hr over broken, rough terrain. Our testing showed moderate success in our suspension design. We successfully damped the force being transmitted to both the front and rear quarter car somewhat, while ensuring that the modes of vibration fell into the desired shapes at the desired frequency ranges.
\end{abstract}
\section{Introduction}
\PARstart{T}{he} suspension system of a car heavily affects the stability of the ride, and the comfort of the
driver and passengers. Minimizing impact to cargo is often a concern as well. Not only the magnitude, but the
mode of vibration experienced must be considered. For a comfortable ride, it is generally
desired to have the bounce motion occur at lower frequencies, and the pitching motion occur at
higher frequencies for better handling and control. Of course, the magnitude of these vibrations should be
minimized overall to the greatest extent possible. We used a half-car model for our analysis.
The fixed parameters of the half-car model are given as follows:
\begin{table}[!hbt]
\begin{center}
\caption{Simulation Parameters}
\label{tab:simParameters}
\begin{tabular}{|c|c|}
\hline
Half-car mass & $m=1000$ $ kg$ \\
\hline
Moment of Inertia about CoM & $J_0=240$ \(kg/m^2\) \\
\hline
Wheel Base & $l_{total}=2 m$\\
\hline
\end{tabular}
\end{center}
\end{table}
Further design parameters specify that the vibration of the vehicle be primarily bouncing for frequencies between
1.4 Hz and 2.4 Hz (8.8 to 15.1 rad/sec), and predominantly pitching at frequencies from 2.5 to 3.3 Hz (15.7 to 21. rad/sec). The damping ratio \(\xi_i\) for each quarter car had to fall in the range of .4 to .6. The placement of the center of mass was left to our discretion. Over the course of our analysis we ended up using different damping ratios for the front and the rear. We also had to change the location of our center of mass from its initial position.
Our dynamic analysis of the problem yielded two equations of motion, from which we created a matrix equation. Algebraic manipulation of the matrix equation allowed us to derive equations for the amplitudes of bouncing and pitching motion. We assumed front and rear damping ratios, and an arbitrary center of mass. Our numerical modeling showed us that our initially chosen values would not be appropriate, and through a trial and error process, we settled on our final values for spring and damping constants, and settled on a center of mass slightly to the rear of the center.
The second phase of our research involved testing our vehicle traveling through conditions simulating broken terrain. We sought to determine the maximum bouncing and pitching displacements experienced by the car. We determined that what we were concerned with was the steady state response of the vehicle, since the long term vibrations as the vehicle is in motion are what we are trying most to limit. Our steady state response was developed from many small calculations, but ultimately allowed us to describe the car's response to harmonic forcing.
\begin{figure}[h!]
\includegraphics[scale=.5]{problem}
\caption{Diagram of half-car.}
\end{figure}
\section{Analysis}
We began our analysis by drawing a free-body diagram for the entire half-car. This
yielded two equations of motion in terms of x and \(\theta\).
\[m\ddot{x} = -k_1(x-l_1\theta)-k_2(x+l_2\theta)-c_1(\dot{x}-l_1\dot{\theta})-c_2(\dot{x}-l_2\dot{\theta})\] % X-Direction EoM
\[J_0\ddot{\theta}=k_1l_1(x-l_1\theta)-k_2l_2(x+l_2\theta)+c_1l_1(\dot{x}-l_1\dot{\theta})-c_2l_2(\dot{x}-l_2\dot{\theta})\] % Angular EoM
\begin{table}[!hbt]
\begin{center}
\caption{Definitions of Variables}
\label{tab:defVariable}
\begin{tabular}{|c|c|}
\hline
Distance from rear, front axles to CoM & $l_1, l_2$ \\
\hline
Spring constant for rear, front axles& $k_1, k_2$ \\
\hline
Damping constant for rear, front axles& $c_1, c_2$\\
\hline
Damping ratio, for rear, front axles& \(\xi_i=b_i/{2m_iw_{ni}} \)\\
\hline
\end{tabular}
\end{center}
\end{table}
Using the following relations,
\begin{center}\(x=\bar{X}e^{st}\)\qquad and\qquad\(\theta=\bar{\Theta}e^{st}\)
\(\dot{x}=s\bar{X}e^{st}\)\qquad and\qquad\(\dot{\theta}=s\bar{\Theta}e^{st}\)
\(\ddot{x}=s^2\bar{X}e^{st}\)\qquad and\qquad\(\ddot{\theta}=s^2\bar{\Theta}e^{st}\)
\end{center}
we were able to put together a matrix expression for \(\bar{X}\) and \(\bar{\Theta}\) in terms of the complex variable \(s=\sigma + j\omega\).
\[\begin{bmatrix}
A & B\\
C & D
\end{bmatrix}
\begin{bmatrix}
\bar{X} \\
\bar{\Theta}
\end{bmatrix}
=\vec{0}
\]
where
\\
\qquad\(A=ms^2+(c_1+c_2)s+(k_1+k_2)\)
\qquad\( B=(c_cl_2-c_1l_1)s+(k_2l_2-k_1l_1)\)
\qquad\( C=(c_2l_2-c_1l_1)s+(k_2-k_1)\)
\qquad\( D=J_0s^2+(c_1l_1^2+c_2l_2^2)s+(k_1l_1^2+k_2l_2^2)\)\\
At this point, we now had an expression to evaluate the bouncing and pitching modes of vibration, \((\frac{\bar{X}}{\bar{\Theta}})\). Using the assumption the real component of \(s\) is negligible, i.e. \(s=j\omega\), we used C++ code to iterate through the full range of frequencies with which we were concerned, from 8.8 to 21 radians/sec. This necessitated that we choose an arbitrary damping ratio for the front and rear that we chose to set at \(\xi_1=\xi_2=.4\). At this point we assumed that the real component of s was negligible. That is to say, we assumed \(s=j\omega\). %We also at this time chose to arbitrarily set the damping ratio for each quarter car at \(\xi_i=.4\).
Initially, we set the center of mass at .8 meters from the front axle arbitrarily. Using C++ code, we solved for spring and damper values for the front and rear that exhibited appropriate bouncing and pitching behavior at the desired frequencies. Using MATLAB, however, we returned to our earlier calculations and solved for the roots of the determinant of the impedance matrix, four values of $s$, in two complex conjugate pairs. We discovered that due to our earlier assumption, \(s=j\omega\), one pair of complex conjugate roots was outside our desired frequency range. We then adjusted our inputs to the C++ code, and using a trial and error process ultimately changed the damping ratio of the front axle to \(\xi_2 = .5\) and set the center of mass 1.2 meters from the front. This yielded spring and damping constants which met our design specifications.
\begin{table}[!hbt]
\begin{center}
\caption{Final chosen physical constants}
\label{tab:kc}
\begin{tabular}{|c|c|}
\hline
Front spring constant & \(k_2 = 250.0 \text { }kN/m\) \\
\hline
Front damping constant & \(c_2=10\text { }kN \cdot s/m\) \\
\hline
Rear spring constant & \(k_1 = 464.64\text { } kN/m\)\\
\hline
Rear damping constant & \(c_1 = 4.224\text { }kN \cdot s/m\) \\
\hline
Location of center of mass & 1.2 meters from front of car (\(l_2 = 1.2 \text{ }m\))\\
\hline
\end{tabular}
\end{center}
\end{table}
Using these constants with the C++ code developed earlier, we were able to develop a data set describing the vibrational modes over the frequency spectrum of interest.
We exported this code to an Excel spreadsheet, included in the appendix, from which we produced the graph shown in Fig.~\ref{fig:vibrationalMode}.
\begin{figure}[h!]
\includegraphics[scale=.52]{modes}
\caption{Modes of vibration across the frequency spectrum}
\label{fig:vibrationalMode}
\end{figure}
The second part of our analysis involved testing the vehicle to determine the maximum amplitude of the bouncing and pitching experienced by the car traveling at a speed of \(\vec{v}=40\) km/hr over rugged terrain. This terrain was modeled by a harmonic displacement of amplitude .1 meters, over a period of 4 meters. For these calculations we assumed a steady state response. This seemed reasonable, since we were mostly concerned with the ongoing vibrations, and not momentary perturbations. The steady state response for a two degree of freedom system grows out of the ss solution for a single degree of freedom system,
\[ x_{ss}=Y_0[TR]sin(\omega t+\phi)\]
where \[[TR]_i=(\frac{\sqrt{1+(2\xi_ir_i)^2}}{\sqrt{(1+r_i^2)^2+(2\xi_ir_i)^2}} \text{ and } r_i = \frac{\omega}{\omega_{ni}}\]
and \(Y_0\) is the amplitude of the harmonic displacement. In a single degree of freedom system, the \(i\) subscripts are not necessary of course, as there is only one damping ratio, natural frequency, etc.
Having a system with two degrees of freedom gives us two steady state solutions in terms of \(x(t) \text{ and } \theta(t)\), each somewhat more complicated than the equation above. The phase shift of \(\pi\) occurs in the second term of each equation as a result of the length of the wheel base being 2 meters to the 4 meter period of the road surface.
\[ x_{ss} = Y_0([TR]_2sin(\omega t +\phi_1+\alpha_1)+[TR]_1sin(\omega t
+ \pi +\phi_2 + \alpha_2))\]
\[\theta_{ss}=Y_0(\frac{-[TR]_2}{l_2}sin(\omega t+\phi_1+\alpha_1)+\frac{[TR]_1}{l_1}sin(\omega t + \pi + \phi_2 + \alpha_2))\]
We used a few standard formulas to find values for \(\phi_1, \phi_2, \alpha_1,\alpha_2, \omega_{n1}, \text{ and } \omega_{n2}\)
\[\phi_i=tan^{-1}(\frac{-2\xi_ir_i}{1-r_1^2}), \text{ } \alpha_i=tan^{-1}(-2\xi_ir_i), \text{ }\omega_{ni}=\sqrt{\frac{k_i}{m_i}} \]
\begin{table}[!ht]
\begin{center}
\caption{Calculated values}
\label{tab:calculated}
\begin{tabular}{|c|c|}
\hline
\(\text{r}_1\) & \(.6981\) \\
\hline
\(\text{r}_2\)& \(1.983\) \\
\hline
\(\text{[TR]}_1\)& \(.7894\)\\
\hline
\(\text{[TR]}_2\)& \(.7423\)\\
\hline
\(\phi_1 \)& \(-.9374 \text{ }rad\)\\
\hline
\(\phi_2\)& \(3.637 \text{ }rad\)\\
\hline
\(\alpha_1\)& \(-.6094 \text{ } rad\)\\
\hline
\(\alpha_2\)& \(-1.008 \text{ }rad\)\\
\hline
\(\omega_{n1}\) & \(25.00\text{ }rad/s\)\\
\hline
\(\omega_{n2}\)& \(8.80 \text{ }rad/s\)\\
\hline
\end{tabular}
\end{center}
\end{table}
Using the information above, we were able to plot the resultant functions in MATLAB and calculate the maximum bouncing and pitching motions in Fig.~\ref{fig:bounceAndPitch}. Our maximum bounce displacement is .0758 meters, or 7.58 cm. The maximum angular displacement is .1407 radians, or about 8.1\(^\circ\).
\begin{figure}[!ht]
\centering
\includegraphics[scale=.55]{extracred}
\caption{Bouncing and pitching motion as a function of time}
\label{fig:bounceAndPitch}
\end{figure}
\section{Conclusion}
The spring and damping constants we arrived at, while large, seem reasonable for a 2000 kilogram vehicle. We were pleased with our efforts to constrain the vibrational motion to the correct modes at the desired frequencies. Given more time to refine our input parameters, we would seek to improve our result by increasing the slope of our curve to achieve even higher mode vibrations(more predominantly bouncing) at the low end of the frequency range, and smaller mode vibrations(more predominantly pitching) at the higher frequencies. Overall, our max bounce mode of \(\frac{\bar{X}}{\bar{\Theta}}=2.02\) and our max pitching mode of \(\frac{\bar{X}}{\bar{\Theta}}=.59\) seem like reasonable results.
While we are not completely thrilled by our car's performance under road conditions, the fact that both transmission ratios are less than zero indicates that our efforts have not been in vain. The maximum pitching motion of 8.1\(^\circ\) seems acceptable for a vehicle traveling over such rough terrain. We were somewhat disappointed with the maximum bounce of 7.58 cm. It seems like we should have been able to reduce the transmitted force more than we succeeded in doing. Despite this, we were otherwise satisfied with our efforts at constraining the vibrational modes to the desired frequency ranges. The vehicle as we designed it seems like a heavily loaded truck, probably for military or industrial use, in light of the rough terrain that we tested it over, and the location of the center of mass behind the geometric center of the truck.
\onecolumn
\textwidth=456pt
\paperwidth=577pt
\hoffset=-30pt
\newpage
\newpage
\clearpage
\pagestyle{headings}
\section{Appendix A}
\centering C++ Code
\footnotesize
\begin{lstlisting}
/*
* File: vibproj.cpp
* Author: Paul
*
* Created on November 12, 2015, 11:15 AM
*/
#include <cstdlib>
#include <stdio.h>
#include <string.h>
#include <math.h>
using namespace std;
/*
*
*/
int main(int argc, char** argv) {
double totalMass = 1000.0;
double xiFront = 0.5;
double xiRear = 0.4;
double length = 2.0;
double frontRatio = .4;
double mFront = totalMass * frontRatio;
double lFront = length - (frontRatio * length);
double wnFront = 8.8;
double kFront = wnFront * wnFront * mFront;
double cFront = 2 * xiFront * wnFront * mFront;
double mRear = totalMass * (1 - frontRatio);
double wnRear = 8.8;
double kRear = wnRear * wnRear*mRear;
double cRear = 2 * xiRear * wnRear*mRear;
double lRear = length - lFront;
double J0 = 240.0;
double wMode = 8.80;
double modeMag = 0;
double max = 0;
double min = 1;
double step = .05;
// printf("%0.2f, %0.2f, %0.2f, %0.2f, %0.2f\n", mFront, lFront, wnFront, kFront, cFront);
while (wMode < 21) {
while (wnRear < 20) {
wnFront = 8.80;
kFront = wnFront * wnFront*mFront;
cFront = 2 * xiFront * wnFront*mFront;
kRear = wnRear * wnRear*mRear;
cRear = 2 * xiRear * wnRear*mRear;
while (wnFront <= 25.05) {
double A = (cRear * lRear - cFront * lFront) * wMode;
double B = kRear * lRear - kFront*lFront;
double E = wMode * wMode * totalMass;
double F = wMode * (cRear + cFront);
modeMag = sqrt(A * A + B * B) / sqrt(E * E + F * F);
if (modeMag > max) {
// printf("wMode: %0.2f\n", wMode);
// printf("Front: Wn: %0.3f\tc: %0.3f\tk: %0.3f\n", wnFront, cFront, kFront);
// printf("Rear: Wn: %0.3f\tc: %0.3f\tk: %0.3f\n", wnRear, cRear, kRear);
// printf("maxMag = %0.3f\n", max);
max = modeMag;
}
if (modeMag < min) {
// printf("wMode: %0.2f\n", wMode);
// printf("Front: Wn: %0.3f\tc: %0.3f\tk: %0.3f\n", wnFront, cFront, kFront);
// printf("Rear: Wn: %0.3f\tc: %0.3f\tk: %0.3f\n", wnRear, cRear, kRear);
// printf("minMag = %0.3f\n", min);
min = modeMag;
}
if (wnFront >= 24.96 && wnRear <= 8.81) {
printf("%0.2f,%0.2f", wMode, modeMag);
printf(",%0.2f,%0.2f,%0.2f,%0.2f\n", kFront, cFront, kRear, cRear);
}
wnFront += step;
kFront = wnFront * wnFront*mFront;
cFront = 2 * xiFront * wnFront*mFront;
}
wnRear += step;
}
wMode += step;
wnFront = 8.80;
kFront = wnFront * wnFront*mFront;
cFront = 2 * xiFront * wnFront*mFront;
wnRear = 8.80;
kRear = wnRear * wnRear*mRear;
cRear = 2 * xiRear * wnRear*mRear;
}
printf(" ================================ ");
return 0;
}
\end{lstlisting}
\normalsize
\pagestyle{headings}
\section{Appendix B}
\centering MATLAB Code
\footnotesize
\begin{lstlisting}
% Part 1
syms m s c1 c2 k1 k2 J l1 l2 t;
clc
m =1000;
l1 = 1.2;
l2 = .8;
J = 240;
c1 = 10000;
k1 = 250000.0;
k2 = 46464;
c2 = 4224;
%s = -37.884 + 18.6275*1i;
A = m.*s.^2 + (c1+c2).*s+(k1+k2);
B = (c2.*l2-c1.*l1).*s + (k2.*l2 - k1.*l1);
C = (c2.*l2-c1.*l1).*s + (k2.*l2 - k1.*l1);
D = J.*s.^2 + (c1.*l1^2 + c2.*l2.^2).*s + (k1.*l1.^2 + k2.*l2.^2);
M = [A B;C D];
%w = 17.4533;
%F = [-.1.*k1.*sin(w.*t) - .1.*k2.*cos(w.*t); .1.*k1.*sin(w.*t)-.1.*k2.*cos(w.*t)];
%inv(M)*F
D = det(M);
P =sym2poly(det(M));
r = roots(P)
A = m.*r.^2 + (c1+c2).*r+(k1+k2);
B = (c2.*l2-c1.*l1).*r + (k2.*l2 - k1.*l1);
mag = abs(-B/A)
%%
% Part 2
t = 0: .01: 1;
w = 17.4533;
x = .07423.*sin(w.*t-1.547)+.07894.*sin(w.*t-3.654);
theta = -.06186.*sin(w.*t-1.547)+.09868.*sin(w.*t-3.654);
xmax = max(x)
thetamax = max(theta)
hold on
xlabel('Time (s)')
plot (t,x,'b')
plot(t, theta,'r'), legend('x(t) in meters', '\theta(t) in radians')
\end{lstlisting}
\normalsize
\pagestyle{headings}
\section{Appendix C}
\centering Specific Frequency Data
\csvautolongtable{vibProj.csv}
\end{document}
