A Simple VBA Tutorial: Frequency response of a mass spring damper

1) The Theory :-

mx”(t)+cx'(t)+kx(t)=f(t)

Taking lapace transforms to obtain the transfer function in complex plain:-

ms^2X(s)+csX(s)+kX(s)=F(s)

assuming 0 initial conditions.

The transfer function is the ratio of output to input (X(s)/F(s) giving us :-

G(s) = X(s)/F(s) = 1/(ms^2+cs+k)

to map this to the frequency domain we replace s with jw (j = -1 complex number)

G(jw) = 1/(-mw^2+cjw+k)

G(jw) = 1/(cjw-mw^2+k)

to split this into real and imaginary part we multiply by the complex conjugate

G(jw) = 1/(cjw-(mw^2-k)) * (cjw+(mw^2-k))/(cjw+(mw^2-k))

G(jw) = (cjw+(mw^2-k)/(-c^2w^2 – (mw^2-k)^2)

splitting this gives us the real and imaginary parts

G(jw) = cjw/(-c^2w^2 – (mw^2-k)^2)    +     (mw^2-k)/(-c^2w^2 – (mw^2-k)^2)

The magnitude of the vibration is the magnitude of G(jw)

The following program calculates the values of G(jw) through the frequency range and plots the magnitude and phase in excel.

2) The code :-

Option Explicit

‘Globals

Dim M As Double

Dim C As Double

Dim k As Double

‘*************************************************************************

‘Program to calculate the response of a mass spring damper to

‘sin input

‘*************************************************************************

Private Sub CommandButton1_Click()

Dim I As Double

Dim R As Double

Dim W As Double

Dim CW As Double

Dim MWK As Double

Dim iW As Integer

Dim Mag As Double

Dim Mag0 As Double

Dim iRow As Integer

Dim Phi As Double

iRow = 2

‘*************************************************************************

‘ Problem setup

‘*************************************************************************

W = 0

‘Read mass damping and stiffness from worksheet

M = Worksheets(1).Cells(2, 6)

C = Worksheets(1).Cells(3, 6)

k = Worksheets(1).Cells(4, 6)

‘Write headings for output in excel

Worksheets(1).Cells(1, 1) = “W”

Worksheets(1).Cells(1, 2) = “Mag/Mag0”

‘*************************************************************************

‘ Calculation

‘ loop 1000 times

‘*************************************************************************

For iW = 0 To 1000

CW = C * W

MWK = M * W * W – k

‘Real and inaginary part of the FRF (Fequency response function)

I = CW / (-CW ^ 2 – MWK ^ 2)

R = MWK / (-CW ^ 2 – MWK ^ 2)

‘Magnitude of FRF

Mag = (I * I + R * R) ^ 0.5

If iW = 0 Then

‘initial magnitude at zero rotational frequency (static response)

Mag0 = Mag

End If

Worksheets(1).Cells(iRow, 1) = W

‘Magifcation ratio

Worksheets(1).Cells(iRow, 2) = Mag / Mag0

‘Phase

Worksheets(1).Cells(iRow, 3) = 180 * (WorksheetFunction.Atan2(R, I)) / 3.145

‘increment W by 0.01

W = W + 0.01

iRow = iRow + 1

Next iW

End Sub

Home ] [ Up ]