Message Boards Message Boards

Dsolve to find out the resonance frequency of mass-spring-damper system?

Posted 6 years ago

Hi, I want to use Dsolve to find out the system resonance frequency expressed by systems parameters. For example, for a simple mass-spring-damper system, the solution using Dsolve contains a term like this enter image description here from which I know the resonance frequency is enter image description here.

Now I want to find out the resonance frequency formula of the following system.

enter image description here

F=s*t, s is constant. The system equations can be expressed as differential equations shown below. enter image description here

DSolve[{m1* x''[t] + c1*x'[t] + ks*x[t] + klc*(x[t] - y[t]) == 0, m2 *y''[t] + c2 *y'[t] + klc*( y[t] - x[t]) + ksb* y[t] == s*t}, {x[t], y[t]}, t]

The result I got is very long and strange. It contains terms like enter image description here I don't know what the # means.

And also in the solution there is a t1 term, I don't know where the t1 comes from.

Anyone's any suggestion will be appreciated.

POSTED BY: Chengjun Li
4 Replies

Chengjun,

The #1 terms are likely inside of a RootSum or Root function. They are arguments of a pure function that would need to be solved to get your answer. Mathematica solves your differential equation in terms of the roots of a polynomial expression. If you put numbers in for the parameters you can get the numerical answer. You can also use ToRadicals[] to get the expression in radical form but in your case it will be 4th order and quite messy (which is why MMA defaults to Root for higher order roots).

Regards,

Neil

POSTED BY: Neil Singer

Chengjun,

If your goal is to only get the eigenvalues (which gives you the damped natural frequencies), I would suggest that you create a state space model and find the eigenvalues of the A matrix. In MMA:

In[27]:= ssm = 
  StateSpaceModel[{m1*x''[t] + c1*x'[t] + ks*x[t] + 
      klc*(x[t] - y[t]) == 0, 
    m2*y''[t] + c2*y'[t] + klc*(y[t] - x[t]) + ksb*y[t] == 
     ff[t]}, {{x[t], 0}, {y[t], 0}}, {ff[t]}, {x[t], y[t]}, t];

In[26]:= {a, b, c, d} = Normal[ssm]

Out[26]= {{{0, 1, 0, 0}, {-((klc + ks)/m1), -(c1/m1), klc/m1, 0}, {0, 
   0, 0, 1}, {klc/m2, 
   0, -((klc + ksb)/m2), -(c2/m2)}}, {{0}, {0}, {0}, {1/m2}}, {{1, 0, 
   0, 0}, {0, 0, 1, 0}}, {{0}, {0}}}

In[19]:= eig = Eigenvalues[a]

Out[19]= {(1/(m1 m2^2))
 Root[klc ks m1^3 m2^7 + klc ksb m1^3 m2^7 + 
    ks ksb m1^3 m2^7 + (c1 klc m1^2 m2^5 + c2 klc m1^2 m2^5 + 
       c2 ks m1^2 m2^5 + c1 ksb m1^2 m2^5) #1 + (c1 c2 m1 m2^3 + 
       klc m1^2 m2^3 + ksb m1^2 m2^3 + klc m1 m2^4 + 
       ks m1 m2^4) #1^2 + (c2 m1 m2 + c1 m2^2) #1^3 + #1^4 &, 1], (1/(
 m1 m2^2))Root[
  klc ks m1^3 m2^7 + klc ksb m1^3 m2^7 + 
    ks ksb m1^3 m2^7 + (c1 klc m1^2 m2^5 + c2 klc m1^2 m2^5 + 
       c2 ks m1^2 m2^5 + c1 ksb m1^2 m2^5) #1 + (c1 c2 m1 m2^3 + 
       klc m1^2 m2^3 + ksb m1^2 m2^3 + klc m1 m2^4 + 
       ks m1 m2^4) #1^2 + (c2 m1 m2 + c1 m2^2) #1^3 + #1^4 &, 2], (1/(
 m1 m2^2))Root[
  klc ks m1^3 m2^7 + klc ksb m1^3 m2^7 + 
    ks ksb m1^3 m2^7 + (c1 klc m1^2 m2^5 + c2 klc m1^2 m2^5 + 
       c2 ks m1^2 m2^5 + c1 ksb m1^2 m2^5) #1 + (c1 c2 m1 m2^3 + 
       klc m1^2 m2^3 + ksb m1^2 m2^3 + klc m1 m2^4 + 
       ks m1 m2^4) #1^2 + (c2 m1 m2 + c1 m2^2) #1^3 + #1^4 &, 3], (1/(
 m1 m2^2))Root[
  klc ks m1^3 m2^7 + klc ksb m1^3 m2^7 + 
    ks ksb m1^3 m2^7 + (c1 klc m1^2 m2^5 + c2 klc m1^2 m2^5 + 
       c2 ks m1^2 m2^5 + c1 ksb m1^2 m2^5) #1 + (c1 c2 m1 m2^3 + 
       klc m1^2 m2^3 + ksb m1^2 m2^3 + klc m1 m2^4 + 
       ks m1 m2^4) #1^2 + (c2 m1 m2 + c1 m2^2) #1^3 + #1^4 &, 4]}

Note that this gives you the four eigenvalues which correspond to the four roots of a polynomial (two complex pairs of poles). You can get the radical form with ToRadicals[] but I think it may be more useful to put numerical values in at this point and get frequencies as the radical equations are very long even if simplified.

Regards,

Neil

POSTED BY: Neil Singer
Posted 6 years ago

Dr. Neil,

Thanks a lot for your reply. I really appreciate it. I have got the numerical values of the system resonance frequencies before from the bode plot of the system. Now I have another way to calculate them with your help. It would be better if I can find out the symbolic expression of the resonance frequency so that I can figure out which parameter will influence the frequency and in what manner. But I realize that it might be very difficult to do that now. Thanks again for your kind help!

Chengjun

POSTED BY: Chengjun Li

Chengjun,

In reference to your comment

It would be better if I can find out the symbolic expression of the resonance frequency so that I can figure out which parameter will influence the frequency and in what manner.

You can do this by realizing that the origin of the equations is a mecahnical vibratory system. With this knowledge, you should set the c1 and c2 variables to zero because their only affect is to add damping to the system. You can find the undamped natural frequency of the system and see the affect of your variables. Let's look at only the first eigenvalue -- they are the same except for a sign change:

In[7]:= eig[[1]]

Out[7]= (1/(m1 m2^2))Root[
 klc ks m1^3 m2^7 + klc ksb m1^3 m2^7 + 
   ks ksb m1^3 m2^7 + (c1 klc m1^2 m2^5 + c2 klc m1^2 m2^5 + 
      c2 ks m1^2 m2^5 + c1 ksb m1^2 m2^5) #1 + (c1 c2 m1 m2^3 + 
      klc m1^2 m2^3 + ksb m1^2 m2^3 + klc m1 m2^4 + 
      ks m1 m2^4) #1^2 + (c2 m1 m2 + c1 m2^2) #1^3 + #1^4 &, 1]

In[15]:= natfreq = Simplify[ToRadicals[eig[[1]] /. c1 -> 0 /. c2 -> 0]]

Out[15]= (1/(Sqrt[2] m1 m2^2))(\[Sqrt](-ksb m1^2 m2^3 - ks m1 m2^4 - 
    klc m1 m2^3 (m1 + m2) - Sqrt[
    m1^2 m2^6 (klc^2 (m1 + m2)^2 + 
       2 klc (m1 - m2) (ksb m1 - ks m2) + (ksb m1 - ks m2)^2)]))

In[17]:= Numerator[natfreq^2]

Out[17]= -ksb m1^2 m2^3 - ks m1 m2^4 - klc m1 m2^3 (m1 + m2) - Sqrt[
 m1^2 m2^6 (klc^2 (m1 + m2)^2 + 
    2 klc (m1 - m2) (ksb m1 - ks m2) + (ksb m1 - ks m2)^2)]

In[18]:= Denominator[natfreq^2]

Out[18]= 2 m1^2 m2^4

We know the natural frequency is Sqrt[keff/meff] where keff is the "effective K" of that frequency and meff is the "effective mass" of that frequency. I squared the natural frequency so the keff is the numerator and meff is the denominator. You can use these expressions to see how your parameters affect the K and the M of each frequency and the overall ratio of K/M.

Regards

Neil

POSTED BY: Neil Singer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract