One-Dimensional Modeling and Analysis of Porpoising of Formula 1 Cars using MATLAB

(thumbnail)

- porpoising theory (It includes some of my own guesses so... it could be wrong.)

: casual relationship between y(ride height) and forces from the suspension unit(F_sus) and the aero package(F_aero). The blue dotted line shows the effect of the damper. Note that the directions of F_sus and F_aero are opposite and the work done by each of them has different sign. Also, the spring force(black solid line on the left) is actually non-linear due to the suspension geometry(this has been factored into the modeling.). 

Now we can guess the conditions for porpoising. 

1. Due to some circumstances, floor stall has to be provoked. Inappropriate spring rate and pre-load, bumps on the road, and rolling of the car can be the cause of it.
2. If there is no F_sus that can practically balance F_aero, porpoising must happens. It means the range of F_sus is in between the non-stalled and fully stalled region of F_aero. This always guarantee the underfloor entering stall. It is definately not desirable and therefore the spring rate and pre-load have to be adjusted. Using extreme rocker arm geometry, it might be possible to balance F_sus and stalled F_aero, but it is suspected to require really extreme spring rate, rocker arm radius and pre-load. Most of all, it is never a good idea to keep the underfloor stalled.
3. Even if F_sus and F_aero can be matched, porpoising still can be triggered. If there is any chance for the absolute value of the work done by damper and the work done by aero package, which are the resistance and the power source of porpoising, to be matched through a cycle, porpoising can happen. If the work(absolute value) done by damper always exceed the work done by the aero package, it will be damped through some cycles. 

- about simulation

: polynomials used to model the downforce characteristics when the floor is stalled and not stalled.

One may argue that porpoising should be modeled at least in two-dimension considering front suspension movement, but here is my opinion.

1. Observing actual footage, whole car almost pivots around front axle while porpoising and it seems that there is not much front suspension movement.
2. The frequency of porpoising is quite consistent and no weird interaction between front and rear suspension was observed. 

- simulation results

: stall starting-ending: 0.06-0.08 m, intended ride height: 0.055 m, kdmp=2000 N/(mps)

Above is the condition 2. mentioned in the beginning of this article. 

: stall starting-ending: 0.06-0.08 m, intended ride height: 0.065 m, kdmp=4200 N/(mps), porpoising triggered(1 bump, 0.02 m, 10 Hz(Note that speed is 85m/s)) but damped through cycles

Above is the condition 3.

: stall starting-ending: 0.06-0.08 m, intended ride height: 0.065 m, kdmp=2000 N/(mps), porpoising triggered(1 bump, 0.02 m, 10 Hz)

Only kdmp is different compared to the condition 3. but porpoising is continuously happening. Note that it can endure some small disturbance not causing floor stall like the very beginning. 

The condition 3. and this may imply the track dependency of porpoising(e.g. Mercedes W13). There are two possibilities.

1. Porpoising can be damped through some cycles, but many bumps keep triggering it.
2. Porpoising is triggered by a bump and does not damped enough.

- effects of aerodynamic characteristics, spring rate, spring pre-load(ride height), and damping rate

case	Base	1. Relaxed Stall Hysterisis	2. Decreased  Spring Rate	3. Increased Ride Height	4. Increased Damping Rate
displacement frequency (1/s)	4.03	4.31	3.64	3.80	3.96
amplitude total (m)	0.034873	0.023384	0.036592	0.035805	0.02677
amplitude upper (m)	0.091043	0.081983	0.09335	0.094815	0.086191
amplitude lower (m)	0.05617	0.058600	0.056758	0.05901	0.059421
acceleration amplitude total  (m/s^2)	24.46	19.93	22.01	24.40	23.65
acceleration amplitude upper (m/s^2)	13.45	12.39	12.70	14.49	12.56
acceleration amplitude lower (m/s^2)	-11.01	-7.54	-9.31	-9.91	-11.09

case_Base: (stall starting-ending: 0.06-0.08 m, kspr=800000 N/m, rdh=0.065 m, kdmp=2000 N/(mps), porpoising triggered and continued)

Red for increased, blue for decreased factors.

Case 1. is a modeling of reducing stall hysteresis(0.06-0.07 m). Comparing this with the base case, we can observe that both displacement and acceleration amplitude have reduced and the frequency has increased. This agrees well with the F1TV Tech Talks data. 

Reducing the amplitude of displacement and acceleration is indeed very desirable and we can understand why every team is so obsessed with floor design and cut-outs. In addition of my personal note, I believe the rounded Redbull underfloor and lots of corotating vortex generating shape can also help reducing stall hysteresis.

If we completely remove stall hysteresis(e.g. 0.06-0.06 m), porpoising is removed or damped, which is understandable.

The spring rate is decreased in the case 2. As a result, the displacement amplitude has slightly increased and the acceleration amplitude and the frequency have decreased. Note that the minimum and maximum displacement has moved higher.

If we increase the spring rate up to 1200000, exactly opposite happens. 

I should mention that spring rate is not that free to adjust as it is related to many other factors such as vehicles attitude in acceleration, braking, bumps and speed-ride height relation. 

By adjusting the spring pre-load only, ride height is increased up to 0.07 m. We can see that the acceleration amplitude and frequency has reduced while the displacement amplitude increasing slightly. So this can be another way to improve porpoising. 

If we raise the ride height up to 0.075 m, the porpoising is not triggered with the bump.

Of course, raising ride height is not good for overall aerodynamics.

Increasing damping rate can improve every parameters, obviously. However, I should mention that damping rate is closely related to vehicle dynamics, which means there will be some compromises.

[MATLAB code]

Aware that I'm just a university student now and there can be lots of inaccurate coefficients and errors below.

[representative code - porpoising triggered by a bump, ride height = 0.065 m]

! All rights reserved !

% porpoising simulator

% numerics

delt=0.0001;

yini=0.07;

vini=0;

tend=10;

exmmod=1;

texms=6;

std=4;

exms=int64(fix(texms/delt));

itr=int64(fix(tend/delt));

bmpmod=1;

bmps=(itr/tend)*(2);

bmpe=(itr/tend)*(2+0.05);

Ab=0.02;

k=((2*pi)/(double(itr)/tend))*(10);

kela=0.8; % COR

% suspension geometry

th1=(pi/180)*(-35);

th2=(pi/180)*(0);

th3=(pi/180)*(35);

r1=0.08;

r2=0.1;

r3=0.6;

kspr=800000;

kdmp=2000;

rdh=0.065;

% C.G. info

m=795;

g=9.8;

kdst=(65/100);

% aero info

p=1.225;

Cl=4.5;

A=2*0.5;

V=85;

a0=(-0.452);

b0=(47.4);

c0=(-539);

d0=(2520);

e0=(-4210);

as=(0.6);

bs=(0.0274);

cs=(9.11);

ds=(-202);

es=(1790);

stalls=0.06;

stalle=0.08;

rdhref=0.1;

ddec=(0.1/rdhref);

% spring preload

Fdwn=(0.5*p*Cl*A*(V^2))*(1-ddec*(rdh-rdhref))*(a0+b0*rdh+c0*(rdh^2)+d0*(rdh^3)+e0*(rdh^4));

%Fdwn=0;

sprpre=(kdst*(m*g+Fdwn)*r1*sin(pi/2-th3-th1))/(2*kspr*r2*cos(th2)*sin(th3));

% -------     END OF SET-UP SECTION     -------

t=linspace(0,tend,itr);

y=zeros(1,itr);

v=zeros(1,itr);

a=zeros(1,itr);

ydiff=zeros(1,itr);

Fw=zeros(1,itr);

Fd=zeros(1,itr);

Fspr=zeros(1,itr);

Fdmp=zeros(1,itr);

Fnet=zeros(1,itr);

y(1)=yini;

v(1)=vini;

a(1)=0;

% before run disp

fprintf('spring preload: %f m \n',sprpre);

fprintf('iteration scheduled: %f \n',itr);

cnt=1;

stl=0;

bttmitr=0;

for i=2:itr

    

    Fw(i-1)=kdst*m*g;

    Fdtemp=(0.5*p*Cl*A*(V^2))*kdst;

    if y(i-1)<0.2

        if (y(i-1)<stalls)||((stl==1)&&(y(i-1)<stalle))

            Fd(i-1)=Fdtemp*(1-ddec*(y(i-1)-rdhref))*(as+bs*y(i-1)+cs*(y(i-1)^2)+ds*(y(i-1)^3)+es*(y(i-1)^4));

            stl=1;

        else

            Fd(i-1)=Fdtemp*(1-ddec*(y(i-1)-rdhref))*(a0+b0*y(i-1)+c0*(y(i-1)^2)+d0*(y(i-1)^3)+e0*(y(i-1)^4));

            stl=0;

            

        end

    else

        Fd(i-1)=Fdtemp*(1-ddec*(0.2-rdhref))*0.91;

    end

    Fspr(i-1)=2*kspr*(r2*(sin(th2+(rdh-y(i-1))/r3)-sin(th2))+sprpre)*r2*cos(th2+(rdh-y(i-1))/r3)*sin(th3)/(r1*sin(pi/2-th3-th1-(rdh-y(i-1))/3));

    Fdmp(i-1)=2*kdmp*v(i-1);

    Fnet(i-1)=Fspr(i-1)-Fw(i-1)-Fd(i-1)-Fdmp(i-1);

    a(i)=Fnet(i-1)/(m*kdst);

    v(i)=v(i-1)+a(i)*delt;

    y(i)=y(i-1)+v(i-1)*delt+0.5*a(i)*(delt^2);

    

    % bump sim

    if (bmpmod==1)&&((i>=bmps)&&(i<=bmpe))

        ydiff(i)=Ab*sin(k*double(i));

    end

    y(i)=y(i)-(ydiff(i)-ydiff(i-1));

  

    % bottoming sim

    if (y(i)<0)

        %{

        if bttmitr==0

            disp(i)

            disp(ydiff(i))

            disp(y(i))

            disp(v(i))

            disp(v(i)-(ydiff(i)-ydiff(i-1))/delt)

        end

        %}

        

        temp=v(i)-(ydiff(i)-ydiff(i-1))/delt;

        v(i)=(-1)*temp*kela;

        a(i)=(v(i)-temp)/delt;

        y(i)=0;

        

        bttmitr=bttmitr+1;

    end

    

    cnt=cnt+1;

end

% post processing

pycnt=0;

if exmmod==1

    mod=0;

    chku=0;

    chkl=0;

    memtu=zeros(1,std);

    memtl=zeros(1,std);

    memsu=zeros(1,std);

    memsl=zeros(1,std);

    prdsum=0;

    

    % displacement examination

    for i=exms:itr

        if y(i-1)<y(i)

            switch mod

                case 0

                    mod=1;

                case (-1)

                    mod=1;

                    chkl=chkl+1;

                    memtl(chkl)=(i-1);

                    memsl(chkl)=y(i-1);

                otherwise

                    mod=1;

            end

        else

            switch mod

                case 0

                    mod=(-1);

                case 1

                    mod=(-1);

                    chku=chku+1;

                    memtu(chku)=(i-1);

                    memsu(chku)=y(i-1);

                otherwise

                    mod=(-1);

            end

        end

    

    

        if (chku==std)&&(chkl==std)

            break;

        end

        pycnt=pycnt+1;

    end

    if (chku~=std)||(chkl~=std)

        fprintf('not enough local maxima and local minima of y to meet the examination standard.');

    end

    for j=2:std

        prdsum=prdsum+delt*((memtu(j)-memtu(j-1))+(memtl(j)-memtl(j-1)));

    end

    

    frqs=(2*(std-1))/prdsum;

    ampu=mean(memsu);

    ampl=mean(memsl);

    

    

    % resetting variables

    mod=0;

    chku=0;

    chkl=0;

    memtu=zeros(1,std);

    memtl=zeros(1,std);

    memau=zeros(1,std);

    memal=zeros(1,std);

    prdsum=0;

    

    % acceleration examination

    pacnt=0;

    for i=exms:itr

        if a(i-1)<a(i)

            switch mod

                case 0

                    mod=1;

                case (-1)

                    mod=1;

                    chkl=chkl+1;

                    memtl(chkl)=(i-1);

                    memal(chkl)=a(i-1);

                otherwise

                    mod=1;

            end

        else

            switch mod

                case 0

                    mod=(-1);

                case 1

                    mod=(-1);

                    chku=chku+1;

                    memtu(chku)=(i-1);

                    memau(chku)=a(i-1);

                otherwise

                    mod=(-1);

            end

        end

    

    

        if (chku==std)&&(chkl==std)

            break;

        end

        pacnt=pacnt+1;

    end

    if (chku~=std)||(chkl~=std)

        fprintf('not enough local maxima and local minima of  to meet the examination standard.');

    end

    for j=2:std

        prdsum=prdsum+delt*((memtu(j)-memtu(j-1))+(memtl(j)-memtl(j-1)));

    end

            

    frqa=(2*(std-1))/prdsum;

    ampaccu=mean(memau);

    ampaccl=mean(memal);

end

% after run disp

fprintf('actual iteration: %d \n',cnt);

fprintf('bottoming iteration: %d \n\n',bttmitr);

if exmmod==1

    fprintf('displacement frequency: %f 1/s \n',frqs);

    fprintf('amplitude upper: %f m \n',ampu);

    fprintf('amplitude lower: %f m \n',ampl);

    fprintf('amplitude total: %f m \n',ampu-ampl);

    fprintf('y counter: %d \n\n',pycnt);

    

    fprintf('acceleration frequency: %f 1/s \n',frqa);

    fprintf('acceleration amplitude upper: %f m/s^2 \n',ampaccu);

    fprintf('acceleration amplitude lower: %f m/s^2 \n',ampaccl);

    fprintf('acceleration amplitude total: %f m/s^2 \n',ampaccu-ampaccl);

    fprintf('a counter: %d \n\n',pacnt);

end

subplot(3,1,1);

plot(t,y);

ylim([0, 0.1]);

xlabel('(sec)');

ylabel('(m)');

title('ride height');

subplot(3,1,2);

hold all

plot(t,Fspr+Fdmp);

plot(t,Fd+Fw);

hold off

ylim([0, inf]);

xlabel('(sec)');

ylabel('(N)');

title('net F_{sus} and F_{ext}');

legend('F_{sus}', 'F_{ext}', 'Location', 'eastoutside');

%{

subplot(2,3,2);

plot(t,Fspr);

ylim([0, inf]);

title('Fspring');

%}

%{

subplot(2,3,3);

plot(t,Fd);

ylim([0, 15000]);

title('Fdownforce');

%}

%{

subplot(2,3,4);

plot(t,ydiff);

ylim([-Ab*1.1, Ab*1.1]);

title('gnd fluc');

%}

%{

subplot(2,3,5);

plot(t,Fdmp);

ylim([-inf, inf]);

title('Fdmp');

%}

%%{

subplot(3,1,3);

plot(t,a);

ylim([-15, 15]);

xlabel('(sec)');

ylabel('(m/s^2)');

title('vert. acceleration');

%%}

%{

subplot(2,3,6);

plot(t,Fnet);

ylim([-inf, inf]);

title('Fnet');

%}

[raw data]

>> case_PwTrgrd_Base (stall starting-ending: 0.06-0.08 m, kspr=800000 N/m, rdh=0.065 m, kdmp=2000 N/(mps), porpoising triggered and continued)

spring preload: 0.015727 m 

displacement frequency: 4.030091 1/s 

amplitude upper: 0.091043 m 

amplitude lower: 0.056170 m 

amplitude total: 0.034873 m 

y counter: 9430 

acceleration frequency: 4.030091 1/s 

acceleration amplitude upper: 13.453712 m/s^2 

acceleration amplitude lower: -11.006198 m/s^2 

acceleration amplitude total: 24.459909 m/s^2 

a counter: 9214 

>> case_PwTrgrd_1. (stall starting-ending: 0.06-0.07 m, kspr=800000 N/m, rdh=0.065 m, kdmp=2000 N/(mps), porpoising triggered and continued)

spring preload: 0.015727 m 

displacement frequency: 4.305705 1/s 

amplitude upper: 0.081983 m 

amplitude lower: 0.058600 m 

amplitude total: 0.023384 m 

y counter: 8718 

acceleration frequency: 4.306014 1/s 

acceleration amplitude upper: 12.393746 m/s^2 

acceleration amplitude lower: -7.540462 m/s^2 

acceleration amplitude total: 19.934208 m/s^2 

a counter: 8566 

>> case_PwTrgrd_2. (stall starting-ending: 0.06-0.08 m, kspr=600000 N/m, rdh=0.065 m, kdmp=2000 N/(mps), porpoising triggered and continued)

spring preload: 0.020969 m 

displacement frequency: 3.639010 1/s 

amplitude upper: 0.093350 m 

amplitude lower: 0.056758 m 

amplitude total: 0.036592 m 

y counter: 10553 

acceleration frequency: 3.639010 1/s 

acceleration amplitude upper: 12.696526 m/s^2 

acceleration amplitude lower: -9.314362 m/s^2 

acceleration amplitude total: 22.010888 m/s^2 

a counter: 10220 

>> case_PwTrgrd_3. (stall starting-ending: 0.06-0.08 m, kspr=800000 N/m, rdh=0.07 m, kdmp=2000 N/(mps), porpoising triggered and continued)

spring preload: 0.015899 m 

displacement frequency: 3.800836 1/s 

amplitude upper: 0.094815 m 

amplitude lower: 0.059010 m 

amplitude total: 0.035805 m 

y counter: 9803 

acceleration frequency: 3.800836 1/s 

acceleration amplitude upper: 14.493527 m/s^2 

acceleration amplitude lower: -9.909946 m/s^2 

acceleration amplitude total: 24.403473 m/s^2 

a counter: 9685 

>> case_PwTrgrd_4. (stall starting-ending: 0.06-0.08 m, kspr=800000 N/m, rdh=0.07 m, kdmp=3500 N/(mps), porpoising triggered and continued)

spring preload: 0.015727 m 

iteration scheduled: 100000.000000 

actual iteration: 100000 

bottoming iteration: 0 

displacement frequency: 3.957262 1/s 

amplitude upper: 0.086191 m 

amplitude lower: 0.059421 m 

amplitude total: 0.026770 m 

y counter: 8966 

acceleration frequency: 3.957262 1/s 

acceleration amplitude upper: 12.561485 m/s^2 

acceleration amplitude lower: -11.092781 m/s^2 

acceleration amplitude total: 23.654267 m/s^2 

a counter: 8608