Bode and Nyquist plots of FRF as function of operational frequency

I have a function for which I need to create a Bode plot and Nyquist Plot. The function describes a rotor, and is as follows. I'v built it up in parts - Hxx is a function of w (omega). k=spring stiffness, c=damping coefficient, l=length of shaft, I=moment of inertia, J=Diametral Inertia, W=rotational speed(rad/s), w=an array of frequencies between w=0:100:30000

n1=k*(l^2);
n2=w.*(i*c*(l^2));
n3=w.^2.*I;

d1=(k^2)*(l^2);
d2=w.*(2*i*c*k*(l^2));
d3=(w.^2).*((c^2*l^2)+(2*I*k)+(W*J/l)^2);
d4=(w.^3).*(2*i*c*I);
d5=(w.^4).*((I/l)^2);

Hxx=(n1+n2-n3)./(d1+d2-d3-d4+d5);

I can happily perform the following plots 1. plot(w, Hxx);

and

realHxx=real(Hxx);
imagHxx=imag(Hxx);
% Calculate the magnitude
mag=sqrt(realHxx.^2+imagHxx.^2);
plot(w, mag);

but I am now stuck as to how to produce a Bode and Nyquist plot. I suspect it is something with how to convert the function into an LTI model using frd, but I'm simply out of my depth with the documentation and cannot see how to proceed. Can anyone help please.

Don Howard

Answer 1

[Arkadiy Turevskiy]

You can simply create an frd object:

sys=frd(Hxx,w);
bode(sys,w);
nyquist(sys,w);

Another option is instead of manually calculating the frequency response like you do, simply create a transfer function describing your system.

s=tf('s'); %laplace transform variable

Now modify your code to replace w with s. Also get rid of dot operations, you do not need them anymore:

n1=k*(l^2);
n2=s*(c*(l^2));
n3=(s^2)*I;
d1=(k^2)*(l^2);
d2=s*(2*c*k*(l^2));
d3=(s^2)*((c^2*l^2)+(2*I*k)+(W*J/l)^2);
d4=(s^3)*(2*c*I);
d5=(s^4)*((I/l)^2);
Hxx=(n1+n2+n3)./(d1+d2+d3+d4+d5);

I am not sure if W in d3 is the same as w. If it is, replace with s as well.

Now Hxx is a transfer function object.

Bode plot:

bode(Hxx,w);

Nyquist plot:

nyquist(Hxx,w);

HTH. Arkadiy

Answer 2

[uopdon]

Thankyou Arkadiy, that was most useful ~ I cannot tell you how much. Excellent.

I must add that creating the FRD object gave me the answer I was looking for, whereas the transfer function method yields different results. I had got as far as this yesterday:

n0=k*(l^2);      %w^0
n1=(i*c*(l^2));  %w^1
n2=I;            %w^2
n3=0;
n4=0;

d0=(k^2)*(l^2);                   %w0
d1=(2*i*c*k*(l^2));               %w1
d2=((c^2*l^2)+(2*I*k)+(W*J/l)^2); %w2
d3=(2*i*c*I);                     %w3
d4=((I/l)^2);                     %w4

num=[-n2 n1 n0]
den=[d4 -d3 -d2 d1 d0]

% plot
figure();
Hs=tf(num,den)
bode(Hs,w)

Which is just another way of creating the transfer function. This code provides exactly the same response as yours. I will continue to test it though - I may have made some errors in my programming.

By the way, capital 'W' is different from miniscule 'w': 'W' is the speed the rotor is rotating at (in rad/s) wheras 'w' is a range of frequenices. When Hxx is plotted against 'w' for a given rotational speed, the maximum magnitude indicates the natural frequencies of the shaft.

To repeat,

sys=frd(Hxx,w);
bode(sys,w);
nyquist(sys,w);

is the business - simple like that. Thankyou very much indeed for you help and attention.

sincerely,

Don Howard

Answer 3

[uopdon]

Dear Arkadiy,

The following function tests four different versions of the solution to the problem. You should just be able to copy, paste and run the code thereby seeing the results for yourself. As I mentioned before, only test 3 gives me the 'correct' answer (corresponds with hand calculations and manual method of calculating magnitude). I would be very interested to understand the reasons behind the different results.

 function [] = FRFv2( k, l, c, I, J, W, w )
  %In order to create Bode and Nyquist plots, the Laplace transform is
  %required - or at least an LTI model is required
  %   Detailed explanation goes here

%%announce the method
blah2=sprintf('function: FRFv2( k, l, c, I, J, W )');
disp(blah2);

%%Function test data values

k = 325000
l = 0.2000
c = 10.3330
I = 7.9533e-04
J = 0.0020

W = 3000    % looking at rotor speeds between 0-3000
                   % when 'W' = 0 there is just 1 natural frequency, 'W'>0 then 2 nat freqs

w=0:10:30000   % range of frequencies in which natural frequencies will be found

%%test 1 ~ Calculate Numerator and Denominator coefficients to make a
% transfer function

n0=k*(l^2);          %w^0
n1=(i*c*(l^2));      %w^1
n2=I;                    %w^2
n3=0;
n4=0;

d0=(k^2)*(l^2);                               %w0
d1=(2*i*c*k*(l^2));                          %w1
d2=((c^2*l^2)+(2*I*k)+(W*J/l)^2);   %w2
d3=(2*i*c*I);                                    %w3
d4=((I/l)^2);                                     %w4

num=[-n2 n1 n0];
den=[d4 -d3 -d2 d1 d0];

% plot
figure();
Hs=tf(num,den);
bode(Hs,w);
title('Bode 1 - TF1');
figure();
title('Nyquist 1 - TF1');
nyquist(Hs, w);

%%Test 2 - another method to make a transfer function

s=tf('s'); %laplace transform variable

n1=k*(l^2);
n2=s*(i*c*(l^2));
n3=s^2*I;
d1=(k^2)*(l^2);
d2=s*(2*i*c*k*(l^2));
d3=(s^2)*((c^2*l^2)+(2*I*k)+(W*J/l)^2);
d4=(s^3)*(2*i*c*I);
d5=(s^4)*((I/l)^2);
Hxx=(n1+n2-n3)/(d1+d2-d3-d4+d5);

figure();
bode(Hxx,w);
title('Bode 2 - TF2');
figure();
nyquist(Hxx,w);
title('Nyquist 2 - TF2');

%%Test 3 - using the frd function

n1=k*(l^2);
n2=w.*(i*c*(l^2));
n3=w.^2.*I;

d1=(k^2)*(l^2);
d2=w.*(2*i*c*k*(l^2));
d3=(w.^2).*((c^2*l^2)+(2*I*k)+(W*J/l)^2);
d4=(w.^3).*(2*i*c*I);
d5=(w.^4).*((I/l)^2);

Hxx=(n1+n2-n3)./(d1+d2-d3-d4+d5);

sys=frd(Hxx,w);

figure();
bode(sys,w);
title('Bode3 - frd1');

figure();
nyquist(sys,w);
title('Nyquist3 - frd1');

%%Test 4 - using the Transfer function

s=tf('s'); %laplace transform variable
n1=k*(l^2);
n2=s.*(i*c*(l^2));
n3=s^2.*I;

d1=(k^2)*(l^2);
d2=s.*(2*i*c*k*(l^2));
d3=(s^2).*((c^2*l^2)+(2*I*k)+(W*J/l)^2);
d4=(s^3).*(2*i*c*I);
d5=(s^4).*((I/l)^2);
Hxx4=(n1+n2-n3)/(d1+d2-d3-d4+d5);

figure();
bode(Hxx4,w);
title('Bode 4 - TF');
figure();
nyquist(Hxx4,w);
title('Nyquist 4 - TF');

end

[Arkadiy Turevskiy]

Ok, let me try to explain. In continuous-time domain, transfer functions are defined as functions of Laplace transform variable, s.

You start by writing your transfer function as a function of s. Your transfer function has a numerator of 2nd order and denominator of the 4th order.

So, sys=(a2*s^2+a1*s+a0)/(b4*s^4+b3*s^3+b2*s^2+b1*s+b0).

Important: in this expression, all the coefficients are the real numbers.

Now you want to calculate the magnitude at different frequencies. You do that by substituting s=i*w, where i=sqrt(-1).

Then you write your transfer function as expression of w, not s.

It becomes:

sys=(-a2*w^2+a1*i*w+a0)/(b4*s^4-b3*i*w^3-b2*w^2+b1*i*w+b0).

The problem in your code (and in my original answer, which I corrected now), is that you are trying to define a function of s using coefficients from the expression for w, and that is why you have complex values in denominator and numerator coefficients.

I hope that clarifies it.

So, to make the code correct you do this:

%% test 1 ~ Calculate Numerator and Denominator coefficients to make a

% transfer function

n0=k*(l^2); %w^0

n1=(c*(l^2)); %w^1

n2=I; %w^2

n3=0;

n4=0;

d0=(k^2)*(l^2); %w0

d1=(2*c*k*(l^2)); %w1

d2=((c^2*l^2)+(2*I*k)+(W*J/l)^2); %w2

d3=(2*c*I); %w3

d4=((I/l)^2); %w4

num=[n2 n1 n0];

den=[d4 d3 d2 d1 d0];

% plot

figure();

Hs=tf(num,den);

bode(Hs,w);

title('Bode 1 - TF1');

figure();

title('Nyquist 1 - TF1');

nyquist(Hs, w);

%% Test 2 - another method to make a transfer function

s=tf('s'); %laplace transform variable

n1=k*(l^2);

n2=s*(c*(l^2));

n3=s^2*I;

d1=(k^2)*(l^2);

d2=s*(2*c*k*(l^2));

d3=(s^2)*((c^2*l^2)+(2*I*k)+(W*J/l)^2);

d4=(s^3)*(2*c*I);

d5=(s^4)*((I/l)^2);

Hxx=(n1+n2+n3)/(d1+d2+d3+d4+d5);

figure();

bode(Hxx,w);

title('Bode 2 - TF2');

figure();

nyquist(Hxx,w);

title('Nyquist 2 - TF2');

%% Test 4 - using the Transfer function

s=tf('s'); %laplace transform variable

n1=k*(l^2);

n2=s*(c*(l^2));

n3=s^2*I;

d1=(k^2)*(l^2);

d2=s*(2*c*k*(l^2));

d3=(s^2)*((c^2*l^2)+(2*I*k)+(W*J/l)^2);

d4=(s^3)*(2*c*I);

d5=(s^4)*((I/l)^2);

Hxx4=(n1+n2+n3)/(d1+d2+d3+d4+d5);

figure();

bode(Hxx4,w);

title('Bode 4 - TF');

figure();

nyquist(Hxx4,w);

title('Nyquist 4 - TF');

Now all 4 options give the same result.