Suspension Design Optimization of a 4WD Sedan using MATLAB and EXCEL

This project aims at optimizing the suspension system design for a 4WD Sedan supposed to use a McPherson Strut in the front and a multilink suspension setup in the rear. The program is in its early stages and, this project details the selection process of baseline springs and anti-roll bars using CAE based on some vehicle level targets set by the Design division.

AIM:
To select springs, anti-roll bars, and set up swing arm geometry to meet the targets for ride frequency, roll and pitch gradient while staying within certain constraints.

OBJECTIVES:

1. To calibrate front and rear spring rates.
2. To optimize the front and rear Rollbar diameters.
3. To optimize the horizontal swing arm length and determining the IC locations w.r.t CP.
4. Calculate the front and rear Anti-characteristics w.r.t Acceleration and breaking pitch.
5. Plotting Rear ride frequency and flat ride plot at 65 mph.

ADDITIONAL INFORMATION AND ASSUMPTIONS:

1. It's assumed that spring rates and roll bar diameters are available only in increments of 0.25.
2. It's assumed that the roll bar is solid (inner diameter = 0).
3. Understeer effect is required so a 5% higher load transfer will be added to the front.
4. A 5% deviation in all targets is allowed.
5. The percentages of anti characteristics cannot be negative or greater than 100.

TYPICAL VEHICLE DYNAMICS CASCADE:

ADOPTED WORKFLOW:

1. Ride frequency:

• Calculate the time lag and target rear ride frequency.
• Calibrate the front and rear spring rates using MATLAB.
• Calculate the front and rear wheel rates.
• Calculate the front and rear ride rates.
• Finally, calculate the front and rear actual ride frequencies making sure they reach the target frequencies with minimal error.

2. Anti rollbar:

• Load transfer calculations
• Roll stiffness calculations
• Anti-roll bar design calculations
• Verification

3. Total load transfer without compensation:

• Calculate total load transfer
• Calculate Dynamic front and rear axle weights.
• Calculate front and rear end deflections.
• Calculate pitch angle and pitch gradient.

4. Position the IC location w.r.t CP in horizontal and vertical perspectives using MATLAB such that the values satisfy both the acceleration and braking pitch gradients with a minimum error.

5. Anti-characteristics during Acceleration:

• Calculate front and rear anti angles.
• Calculate front and rear swing arm angles.
• Calculate front and rear anti-lift and anti-squat %.
• Calculate front and rear end deflection due to acceleration.
• Check the actual Acceleration pitch gradient with the target value and note error.

6. Anti-characteristics during Braking:

• Calculate front and rear anti angles.
• Calculate front and rear swing arm angles.
• Calculate front and rear anti-dive and anti-lift%.
• Calculate front and rear end deflection due to braking.
• Check the actual braking pitch gradient with the target value and note error.

GIVEN DATA:

Input Data:

Target Values:

EXCEL CALCULATIONS: (Please go to the end of the page to view the original document)

1. Derived Inputs:

Deliverables:

1. Front spring rate: K_sf = 198.25 lb/in
2. Rear spring rate: K_sr = 273.25 lb/in
3. Front IC vertical distance from CP: h_icf = 14.9367 in ~= 15 in
4. Front IC horizontal distance from CP: l_icf = 78.0672 in ~= 78 in
5. Rear IC vertical distance from CP: h_icr = 14.9367 in ~= 15 in
6. Rear IC horizontal distance from CP: l_icr = 78.0672 in ~= 78 in

2. Ride frequency calculations:

Deliverables:

1. Actual front natural frequency = 1.350 Hz
2. Actual rear natural frequency = 1.527 Hz
3. Front wheel rate: K_wf = 178.92 lb/in
4. Rear wheel rate: K_wr = 174.88 lb/in
5. Front ride rate: K_rf = 160.58 lb/in
6. Rear ride rate: K_rr = 159.11 lb/in

3. Anti Rollbar Calculations:

Deliverables:

1. Front anti rollbar stiffness: s_arbf = 108.306 lb/in
2. Rear anti rollbar stiffness: s_arbr = 15.308 lb/in
3. Front anti rollbar diameter: d_arbf = 18.843 mm
4. Rear anti rollbar diameter: d_arbr = 11.067 mm

4. Acceleration Pitch Calculations :

Deliverables:

1. Front anti-lift %: alf = 13.916 %
2. Rear anti-squat %: ras = 5.964%

5. Braking Pitch Calculations:

Deliverables:

1. Front anti-dive %: fad = 66.039 %
2. Rear anti-lift %: ral = 44.026 %

6. Error Calculation w.r.t target values:

Looking at the above table we can say that all values are satisfying the target values with minimum error and are under the required ranges.

MATLAB CODE:

1.1 Matlab Code to optimize the spring rates w.r.t the target ride frequencies

close all
clear all
clc

% Inputs calculated from excel sheet
m_csf = 862.5; % lb Front corner sprung mass
m_csr = 667.5; % lb Rear corner sprung mass
MR_f = 0.95; % in/in Front Spring motion ratio
MR_r = 0.8; % in/in Rear Spring motion ratio
s_yf = 1567; % lb/in Front Tyre vertical stiffness
s_yr = 1764; % lb/in Rear tyre vertical stiffness
f_ftar = 1.350; % Hz Front ride frequency (target)
f_rtar = 1.527; % Hz Rear ride frequency (target)

% Initializing
K_sf = 0; % lb/in Front spring rate
K_sr = 0; % lb/in Rear spring rate
f_f = ridefeq(K_sf,MR_f,s_yf,m_csf); % Hz Actual front ride frequency
f_r = ridefeq(K_sr,MR_r,s_yr,m_csr); % Hz Actual front ride frequency

% Optimizing the fromt ride frequency
while f_f ~= f_ftar % While the actual front ride frequency not equal to the target front ride frequency
    if f_f > f_ftar % If the actual front ride frequency is greater than the target front ride frequency
        K_sf = K_sf - 0.25; % Decrement the value of the front spring rate with 0.25
    else % If the actual front ride frequency is less than the target front ride frequency
        K_sf = K_sf + 0.25; % Increment the value of the front spring rate with 0.25
    end
    f_f = ridefeq(K_sf,MR_f,s_yf,m_csf); % Recalculate the value with newer variables
end
fprintf('The Front spring rate is %f lb/in @ %f Hz \n',K_sf,f_f);

% Optimizing the rear ride frequency
while f_r ~= f_rtar 
    if f_r > f_rtar
        K_sr = K_sr - 0.25;
    else
        K_sr = K_sr + 0.25;
    end
    f_r = ridefeq(K_sr,MR_r,s_yr,m_csr);
end
fprintf('The Rear spring rate is %f lb/in @ %f Hz \n',K_sr,f_r);

%% Bump Oscillation profile

t = 0.0856643; % Time lag (Calculated in Excel sheet)
time = linspace(0,3,900); % Simulation time

front_amp = zeros(); rear_amp = zeros(); % Initializing the variables for speed

for i = 1:length(time)
    front_amp(i) = exp(-time(i))*sin(2*pi*f_f*time(i)); 
    rear_amp(i) = exp(-time(i))*sin(2*pi*f_r*(time(i)-t));
end

% Flat ride plot
figure()
plot(time,front_amp)
hold on
plot(time,rear_amp)
grid on
title('Flat ride Profile')
xlabel('Time (s)')
ylabel('Amplitude')
legend('Front','Rear')
hold off

1.2 Function ridefeq() code:

function [f] = ridefeq(K_s,MR,s_y,m_cs)
    K_w = K_s*MR^2; % Calculates the wheel rate
    K_r = (K_w*s_y)/(K_w+s_y); % Calculates the ride rate
    f = ((2*pi)^-1)*sqrt(K_r*175.13*2.205/m_cs); % Calculates the ride frequency
    f = round(f*1000)/1000; % Rounding off to the nearest 1000
end

1.3. Output:

Damped flat ride profile:

2.1. MATLAB Code to optimize the Instant centre location w.r.t CP:

close all
clear all
clc 

% Input values
r_lf = 12.27; % in Front tyre loaded radius
r_lr = 12.27; % in Rear tyre loaded radius
h_cg = 17; % in Center of gravity height
l_wb = 98; % in Wheelbase length
b_t = 58; % in Trackwidth
DR_f = 0.7; % Front drive ratio
fed = 0.048; % m Front-end deflection
red = 0.0485; % m Rear-end deflection
Gforce = 1;
BR_f = 0.6; % Rear drive ratio
G_ap = 2.000;
G_bp = 1.000;

% Initializing variables to be optimized
h_icf = 14.9367; % Front IC vertical dist from CP
l_icf = 78.0606; % Front IC horizontal dist from CP
h_icr = h_icf; % Rear IC vertical distance from CP
l_icr = l_icf; % Rear IC horizontal distance from CP

G_apc = G_apcfunc(r_lf,r_lr,h_cg,l_wb,DR_f,fed,red,Gforce,h_icf,l_icf,h_icr,l_icr); % Calculate the G_apc
error_apc = round((100-(G_apc*100/G_ap))*1000)/1000; % Calculate the error

while error_apc ~= 0.000 % While the error in apc is not equal to 0.000
    while G_apc ~= G_ap % While G_apc not equal to G_ap
        if G_apc > G_ap % If G_apc greater than G_ap, do the following
            h_icf = h_icf + 0.0001;
            h_icr = h_icf;
            l_icf = l_icf - 0.0001;
            l_icr = l_icf;
        end
        if G_apc < G_ap % If G_apc lesser than G_ap, do the following
            h_icf = h_icf - 0.0001;
            h_icr = h_icf;
            l_icf = l_icf + 0.0001;
            l_icr = l_icf;
        end
        G_apc = G_apcfunc(r_lf,r_lr,h_cg,l_wb,DR_f,fed,red,Gforce,h_icf,l_icf,h_icr,l_icr); % Recalculate the G_apc with new variables
    end
    
    G_bpc = G_bpcfunc(h_icf,l_icf,h_icr,l_icr,h_cg,l_wb,fed,red,BR_f,Gforce); % Calculate the G_bpc with values of G_apc variables

    while G_bpc ~= G_bp % If G_bpc not equal to G_bp
        if G_bpc > G_bp % If G_bpc is greater do the following
            l_icf = l_icf - 0.0001;
            l_icr = l_icf;
        end
        if G_bpc < G_bp % If G_bpc is lesser than G_bp, do the following
            l_icf = l_icf + 0.0001;
            l_icr = l_icf;
        end
        G_bpc = G_bpcfunc(h_icf,l_icf,h_icr,l_icr,h_cg,l_wb,fed,red,BR_f,Gforce); % Recalculate G_bpc with new variables
    end
    error_apc = round((100-(G_apc*100/G_ap))*10000)/10000; % Recalculate the error in G_apc due to variable alteration in G_bpc
end % While error_bpc = 0.000 exit the loop

G_apc = G_apcfunc(r_lf,r_lr,h_cg,l_wb,DR_f,fed,red,Gforce,h_icf,l_icf,h_icr,l_icr); % Recalculate the G_apc one final time
G_bpc = G_bpcfunc(h_icf,l_icf,h_icr,l_icr,h_cg,l_wb,fed,red,BR_f,Gforce); % Recalculate the G_bpc one final time

% Measure the errors for reference
error_apc = 100-(G_apc*100/G_ap); 
error_bpc = 100-(G_bpc*100/G_bp);

% Print the outputs
fprintf('The front IC vertical distance from CP is %f inn',h_icf);
fprintf('The rear IC vertical distance from CP is %f inn',h_icr);
fprintf('The front IC horizontal distance from CP is %f inn',l_icf);
fprintf('The rear IC horizontal distance from CP is %f inn',l_icr);

2.2. Function G_apcfunc() code:

function [G_apc] = G_apcfunc(r_lf,r_lr,h_cg,l_wb,DR_f,fed,red,Gforce,h_icf,l_icf,h_icr,l_icr)
    alaf = atand((h_icf-r_lf)/l_icf); % deg Front anti-lift angle
    alar = atand((h_icr-r_lr)/l_icr); % deg Rear anti-squat angle
    aswlaf = atand(h_cg/l_wb); % deg Front swing arm angle for 100% antilift
    aswsar = atand(h_cg/l_wb); % deg Front swing arm angle for 100% antilift
    alf = alaf*DR_f*100/aswlaf; % % Front antilift %
    asr = alar*(1-DR_f)*100/aswsar; % % Front antilift %
    fedac = fed*(100-alf)/100; % m Front end deflection due to acceleration 
    redac = red*(100-asr)/100; % m Rear end deflection due to acceleration
    paac = atand((fedac+redac)/(l_wb/39.37)); % deg Pitch angle due to acceleration compensation
    G_apc = paac/Gforce; % Acceleration pitch gradient with compensation
    G_apc = round(G_apc*10000)/10000;
end

2.3. Function G_bpcfunc() code:

function [G_bpc] = G_bpcfunc(h_icf,l_icf,h_icr,l_icr,h_cg,l_wb,fed,red,BR_f,Gforce)
    aadf = atand(h_icf/l_icf);
    aslr = atand(h_icr/l_icr);
    aswadf = atand(h_cg/l_wb);
    aswalr = atand(h_cg/l_wb);
    fad = aadf*BR_f*100/aswadf;
    ral = aslr*(1-BR_f)*100/aswalr;
    fedbc = fed*(100-fad)/100;
    redbc = red*(100-ral)/100;
    abp = atand((fedbc+redbc)/(l_wb/39.37));
    G_bpc = abp/Gforce;
    G_bpc = round(G_bpc*100000)/100000;
end

2.4. Outputs:

*For code explanation please follow the comments in the code itself!

THE LOGIC BEHIND CALIBRATING THE LOCATION OF IC W.R.T CP:

• The ratio between h_icf and l_icf is calibrated to be 5.266 using excel.
  • This implies that h_icf = 5.266 * l_icf, this ratio is different for different suspension setups
  • Since h_icr and l_icr are the same as h_icf and l_icr respectively, we will treat them like the front.
• Let's start with a random value of h_icf. Now calculate the G_apc and its error wrt the target value.
• While the error is not equal to 0.000 we will enter a loop
  • While G_apc is not equal to G_ap, we will enter another loop
    • If G_apc is greater than G_ap,
      • h_icf is incremented by 0.0001 and l_icf is decremented by 0.0001.
      • h_icr is equal to h_icf and l_icr is equal to l_icf like we already discussed.
    • If G_apc is lesser than G_ap,
      • h_icf is decremented by 0.0001 and l_icf is incremented by 0.0001.
      • h_icr is equal to h_icf and l_icr is equal to l_icf.
    • Recalculate G_apc.
    • If G_apc is equal to G_ap with 0% error, exit the loop.
  • Calculate G_bpc with values optimized with G_apc
  • If G_bpc not equal to G_bp, we will enter a loop
    • If G_bpc > G_bp,
      • l_icf is decremented by 0.0001
      • l_icr = l_icf
    • If G_bpc < G_bp
      • l_icf is incremented by 0.0001
      • l_icr = l_icf
    • Recalculate G_bpc, if G_bpc is equal to G_bp, exit the loop
  • Recalculate G_apc with values changed by G_bpc.
  • Calculate the error in G_apc.
    • If the error is not equal to 0.000, go to the start of the super loop.
    • If the error is equal to 0.000, exit all the loops.
• Now you got the optimized values off h_icf, l_icf, h_icr and l_icr such that the G_apc and G_bpc have the lowest errors!

RESOURCES:

• Link to the excel sheet document: Excel document