Simulink files to download
A two degrees of freedom system consisting of two masses connected by springs and subject to 3 different type of input forces is analyzed and simulated using Simulink
The system that is being analyzed is show in the following diagram
In the above, \(F\left ( t\right ) \) is to be taken as each of the following
It is required to find \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) analytically and then to use Matlab’s Simulink software for the analysis.
The mathematical model of the system is first developed and the equation of motions obtained using Lagrangian formulation then the analytical solution is found by solving the resulting coupled second order differential equations for \(m_{1}\) and \(m_{2}\). Next, a simulink model is developed to implement the differential equations and the output \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) from Simulink is shown and compared to the output from the analytical solution.
The following is the free body diagram of the above system
Assuming positive is downwards and that \(x_{2}>x_{1}\), force-balance equations for \(m_{1}\) results in\[ m_{1}\ddot {x}_{1}=F\left ( t\right ) +k_{1}\left ( x_{2}-x_{1}\right ) \] And force-balance equations for \(m_{2}\) results in\[ m_{2}\ddot {x}_{2}=-k_{1}\left ( x_{2}-x_{1}\right ) -k_{2}x_{2}\] Hence the EQM for the system become\begin {align*} m_{1}\ddot {x}_{1}-k_{1}x_{2}+k_{1}x_{1} & =F\left ( t\right ) \\ m_{2}\ddot {x}_{2}+\left ( k_{1}+k_{2}\right ) x_{2}-k_{1}x_{1} & =0 \end {align*}
Or in matrix form\[\begin {pmatrix} m_{1} & 0\\ 0 & m_{2}\end {pmatrix}\begin {pmatrix} \ddot {x}_{1}\\ \ddot {x}_{2}\end {pmatrix} +\begin {pmatrix} k_{1} & -k_{1}\\ -k_{1} & \left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} =\begin {pmatrix} F\left ( t\right ) \\ 0 \end {pmatrix} \] The above can be written in matrix form as\[ M\mathbf {\ddot {x}}+K\mathbf {x}=\mathbf {F}\] Where \(\mathbf {\ddot {x},x,F\,\ }\)are 2 by 1 vectors and \(M\) and \(K\) are the mass and stiffness matrices. The solution to the above is \begin {equation} \mathbf {x}=\mathbf {x}_{h}+\mathbf {x}_{p} \tag {1} \end {equation}
We start by finding \(\mathbf {x}_{h}\) from the following\[\begin {pmatrix} m_{1} & 0\\ 0 & m_{2}\end {pmatrix}\begin {pmatrix} \ddot {x}_{1}\\ \ddot {x}_{2}\end {pmatrix} +\begin {pmatrix} k_{1} & -k_{1}\\ -k_{1} & \left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} =\begin {pmatrix} 0\\ 0 \end {pmatrix} \] Assuming that \(x_{1}\left ( t\right ) =A_{1}\cos \left ( \omega t+\phi \right ) \) and \(x_{2}\left ( t\right ) =A_{2}\cos \left ( \omega t+\phi \right ) \), hence \(\dot {x}_{1}\left ( t\right ) =-\omega A_{1}\sin \left ( \omega t+\phi \right ) \) and \(\dot {x}_{2}\left ( t\right ) =-\omega A_{2}\sin \left ( \omega t+\phi \right ) \) and \(\ddot {x}_{1}\left ( t\right ) =-\omega ^{2}A_{1}\cos \left ( \omega t+\phi \right ) \) and \(\dot {x}_{2}\left ( t\right ) =-\omega ^{2}A_{2}\cos \left ( \omega t+\phi \right ) \). Substituting the above values in the above system results in\[\begin {pmatrix} m_{1} & 0\\ 0 & m_{2}\end {pmatrix}\begin {pmatrix} -\omega ^{2}A_{1}\cos \left ( \omega t+\phi \right ) \\ -\omega ^{2}A_{2}\cos \left ( \omega t+\phi \right ) \end {pmatrix} +\begin {pmatrix} k_{1} & -k_{1}\\ -k_{1} & \left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} A_{1}\cos \left ( \omega t+\phi \right ) \\ A_{2}\cos \left ( \omega t+\phi \right ) \end {pmatrix} =\begin {pmatrix} 0\\ 0 \end {pmatrix} \] Divide by \(\cos \left ( \omega t+\phi \right ) \) since not zero (else no solution exist) we obtain\[\begin {pmatrix} m_{1} & 0\\ 0 & m_{2}\end {pmatrix}\begin {pmatrix} -\omega ^{2}A_{1}\\ -\omega ^{2}A_{2}\end {pmatrix} +\begin {pmatrix} k_{1} & -k_{1}\\ -k_{1} & \left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} A_{1}\\ A_{2}\end {pmatrix} =\begin {pmatrix} 0\\ 0 \end {pmatrix} \] Rewrite the above as\begin {align} \begin {pmatrix} -\omega ^{2}m_{1}A_{1}\\ -\omega ^{2}m_{2}A_{2}\end {pmatrix} +\begin {pmatrix} k_{1}A_{1}-k_{1}A_{2}\\ -k_{1}A_{1}+\left ( k_{1}+k_{2}\right ) A_{2}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \nonumber \\\begin {pmatrix} -\omega ^{2}m_{1}A_{1}+k_{1}A_{1}-k_{1}A_{2}\\ -\omega ^{2}m_{2}A_{2}-k_{1}A_{1}+\left ( k_{1}+k_{2}\right ) A_{2}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \nonumber \\\begin {pmatrix} \left ( -\omega ^{2}m_{1}+k_{1}\right ) A_{1}-k_{1}A_{2}\\ -\omega ^{2}m_{2}A_{2}+\left ( k_{1}+k_{2}\right ) A_{2}-k_{1}A_{1}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \nonumber \\\begin {pmatrix} -\omega ^{2}m_{1}+k_{1} & -k_{1}\\ -k_{1} & -\omega ^{2}m_{2}+\left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} A_{1}\\ A_{2}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \tag {2} \end {align}
From the last equation above, we see that to obtain a solution we must have \[\begin {vmatrix} -\omega ^{2}m_{1}+k_{1} & -k_{1}\\ -k_{1} & -\omega ^{2}m_{2}+\left ( k_{1}+k_{2}\right ) \end {vmatrix} =0 \] since if we had \(\begin {pmatrix} A_{1}\\ A_{2}\end {pmatrix} =0\) then no solution will exist. Therefore, taking the determinant and setting it to zero results in
\begin {align*} \left ( -\omega ^{2}m_{1}+k_{1}\right ) \left ( -\omega ^{2}m_{2}+\left ( k_{1}+k_{2}\right ) \right ) -k_{1}^{2} & =0\\ \allowbreak k_{1}^{2}+k_{1}k_{2}-\omega ^{2}k_{1}m_{1}-\omega ^{2}k_{1}m_{2}-\omega ^{2}k_{2}m_{1}+\omega ^{4}m_{1}m_{2}-k_{1}^{2} & =0\\ \allowbreak \omega ^{4}m_{1}m_{2}+\omega ^{2}\left ( -k_{1}m_{1}-k_{2}m_{1}-k_{1}m_{2}\right ) +\left ( k_{1}^{2}+k_{1}k_{2}-k_{1}^{2}\right ) & =0 \end {align*}
Let \(\omega ^{4}=\lambda ^{2}\), hence the above becomes\[ \lambda ^{2}m_{1}m_{2}+\lambda \left ( -k_{1}m_{1}-k_{2}m_{1}-k_{1}m_{2}\right ) +\left ( k_{1}^{2}+k_{1}k_{2}-k_{1}^{2}\right ) =0 \] Solving for \(\lambda \) gives\begin {equation} \begin {array} [c]{c}\lambda _{1}=\frac {1}{2m_{1}m_{2}}\left ( k_{1}m_{1}+k_{1}m_{2}+k_{2}m_{1}+\sqrt {k_{1}^{2}m_{1}^{2}+k_{1}^{2}m_{2}^{2}+k_{2}^{2}m_{1}^{2}+2k_{1}k_{2}m_{1}^{2}+2k_{1}^{2}m_{1}m_{2}-2k_{1}k_{2}m_{1}m_{2}}\right ) \\ \lambda _{2}=\frac {1}{2m_{1}m_{2}}\left ( k_{1}m_{1}+k_{1}m_{2}+k_{2}m_{1}-\sqrt {k_{1}^{2}m_{1}^{2}+k_{1}^{2}m_{2}^{2}+k_{2}^{2}m_{1}^{2}+2k_{1}k_{2}m_{1}^{2}+2k_{1}^{2}m_{1}m_{2}-2k_{1}k_{2}m_{1}m_{2}}\right ) \end {array} \tag {3} \end {equation} For each of the above solutions, we obtain a different \(\begin {pmatrix} A_{1}\\ A_{2}\end {pmatrix} \) from equation (2) as follows
For \(\lambda _{1}\), (2) becomes\begin {align*} \begin {pmatrix} -\lambda _{1}m_{1}+k_{1} & -k_{1}\\ -k_{1} & -\lambda _{1}m_{2}+\left ( k_{1}+k_{2}\right ) \end {pmatrix}\begin {pmatrix} A_{1}\\ A_{2}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \\\begin {pmatrix} \left ( -\lambda _{1}m_{1}+k_{1}\right ) A_{1}-k_{1}A_{2}\\ -k_{1}A_{1}-\lambda _{1}m_{2}A_{2}+\left ( k_{1}+k_{2}\right ) A_{2}\end {pmatrix} & =\begin {pmatrix} 0\\ 0 \end {pmatrix} \end {align*}
From the first equation above, we have\begin {align*} \left ( -\lambda _{1}m_{1}+k_{1}\right ) A_{1}-k_{1}A_{2} & =0\\ \frac {-\lambda _{1}m_{1}+k_{1}}{k_{1}} & =\left ( \frac {A_{2}}{A_{1}}\right ) ^{\left ( 1\right ) } \end {align*}
Similarly for \(\lambda _{2}\), \[ \frac {-\lambda _{2}m_{1}+k_{1}}{k_{1}}=\left ( \frac {A_{2}}{A_{1}}\right ) ^{\left ( 2\right ) }\] Let\begin {align} r_{1} & =\frac {-\lambda _{1}m_{1}+k_{1}}{k_{1}}=\left ( \frac {A_{2}}{A_{1}}\right ) ^{\left ( 1\right ) }\tag {4}\\ r_{2} & =\frac {-\lambda _{2}m_{1}+k_{1}}{k_{1}}=\left ( \frac {A_{2}}{A_{1}}\right ) ^{\left ( 2\right ) }\nonumber \end {align}
Hence now \(\mathbf {x}_{h}\) can be written as\[\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} _{h}=\begin {pmatrix} A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \\ A_{2}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{2}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \end {pmatrix} \] But \(A_{2}^{\left ( 1\right ) }=r_{1}A_{1}^{\left ( 1\right ) }\) and \(A_{2}^{\left ( 2\right ) }=r_{2}A_{1}^{\left ( 2\right ) }\), hence the above becomes\begin {equation} \begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} _{h}=\begin {pmatrix} A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \\ r_{1}A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +r_{2}A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \end {pmatrix} \tag {5} \end {equation} Now, given numerical values for \(k_{1},k_{2},m_{1},m_{2}\) we can find \(\omega _{1},\omega _{2}\) from (3) above, and next find \(r_{1},r_{2}\) from (4). Hence (5) contains 4 unknowns, \(A_{1}^{\left ( 1\right ) },A_{1}^{\left ( 2\right ) },\phi _{1},\phi _{2}\) which now can be found from initial conditions (after we find the particular solution) which we will now proceed to do.
There are 3 different \(F\left ( t\right ) \) which we are asked to consider
For each of the above, we find \(\mathbf {x}_{p}\) and then add it to \(\mathbf {x}_{h}\) found above in (5) to obtain (1).
Using the standard response for a unit impulse which for a single degree of freedom system is \(x\left ( t\right ) =\frac {1}{m\omega _{n}}\sin \omega _{n}t\), then we write \(\mathbf {x}_{p}\) as\[ \mathbf {x}_{p}=\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} _{p}=\begin {pmatrix} \frac {1}{m\omega _{1}}\sin \omega _{1}t+\frac {1}{m\omega _{2}}\sin \omega _{2}t\\ 0 \end {pmatrix} \] Hence, the general solution becomes\begin {equation} \begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} =\begin {pmatrix} A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \\ r_{1}A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +r_{2}A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \end {pmatrix} +\begin {pmatrix} \frac {1}{m\omega _{1}}\sin \omega _{1}t+\frac {1}{m\omega _{2}}\sin \omega _{2}t\\ 0 \end {pmatrix} \tag {6} \end {equation}
Since unit step is \(1\) for \(t>0\), then, using convolution we write\begin {align*} x_{p}\left ( t\right ) & =\int _{0}^{t}f\left ( \tau \right ) h\left ( t-\tau \right ) d\tau \\ & =\int _{0}^{t}h\left ( t-\tau \right ) d\tau \\ & =\int _{0}^{t}\frac {1}{m\omega _{n}}\sin \omega _{n}\left ( t-\tau \right ) d\tau \\ & =\frac {1}{m\omega _{n}}\left [ \frac {-\cos \left ( \omega _{n}\left ( t-\tau \right ) \right ) }{-\omega _{n}}\right ] _{0}^{t}\\ & =\frac {1}{m\omega _{n}^{2}}\left [ \cos \left ( \omega _{n}\left ( t-\tau \right ) \right ) \right ] _{0}^{t}\\ & =\frac {1}{m\omega _{n}^{2}}\left [ 1-\cos \left ( \omega _{n}t\right ) \right ] \end {align*}
Then, since now we have 2 natural frequencies, we can write \(\mathbf {x}_{p}\) as\[ \mathbf {x}_{p}=\begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} _{p}=\begin {pmatrix} \frac {1}{m\omega _{1}^{2}}\left [ 1-\cos \left ( \omega _{1}t\right ) \right ] +\frac {1}{m\omega _{2}^{2}}\left [ 1-\cos \left ( \omega _{2}t\right ) \right ] \\ 0 \end {pmatrix} \] Hence, the general solution becomes\begin {equation} \begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} =\begin {pmatrix} A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \\ r_{1}A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +r_{2}A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \end {pmatrix} +\begin {pmatrix} \frac {1}{m\omega _{1}^{2}}\left [ 1-\cos \left ( \omega _{1}t\right ) \right ] +\frac {1}{m\omega _{2}^{2}}\left [ 1-\cos \left ( \omega _{2}t\right ) \right ] \\ 0 \end {pmatrix} \nonumber \end {equation}
In this case, we guess that \(x_{1p}=c_{1}\cos \omega t+c_{2}\sin \omega t\), and since there is no forcing function being applied directly on \(m_{2}\) then \(x_{2p}=0\) hence\[ \mathbf {x}_{p}=\begin {pmatrix} x_{1}\left ( t\right ) \\ x_{2}\left ( t\right ) \end {pmatrix} _{p}=\begin {pmatrix} c_{1}\cos \omega t+c_{2}\sin \omega t\\ 0 \end {pmatrix} \] Then \(\dot {x}_{1p}\left ( t\right ) =-\omega c_{1}\sin \omega t+\omega c_{2}\cos \omega t\) and \(\ddot {x}_{1p}\left ( t\right ) =-\omega ^{2}c_{1}\cos \omega t+\omega ^{2}c_{2}\sin \omega t\) and now we substitute these into the original ODE for \(x_{1}\) which is\begin {align*} m_{1}\ddot {x}_{1}-k_{1}x_{2}+k_{1}x_{1} & =F\left ( t\right ) \\ m_{1}\ddot {x}_{1p}-k_{1}x_{2p}+k_{1}x_{1p} & =\sin \omega t \end {align*}
We obtain the following\begin {align*} m_{1}\left ( -\omega ^{2}c_{1}\cos \omega t+\omega ^{2}c_{2}\sin \omega t\right ) -k_{1}\left ( 0\right ) +k_{1}\left ( c_{1}\cos \omega t+c_{2}\sin \omega t\right ) & =\sin \omega t\\ \cos \omega t\left ( -\omega ^{2}c_{1}m_{1}+k_{1}c_{1}\right ) +\sin \omega t\left ( m_{1}\omega ^{2}c_{2}+k_{1}c_{2}\right ) & =\sin \omega t \end {align*}
Hence by comparing coefficients, we obtain\begin {align*} -\omega ^{2}c_{1}m_{1}+k_{1}c_{1} & =0\\ m_{1}\omega ^{2}c_{2}+k_{1}c_{2} & =1 \end {align*}
or\begin {align*} c_{1}\left ( -\omega ^{2}m_{1}+k_{1}\right ) & =0\\ c_{2}\left ( m_{1}\omega ^{2}+k_{1}\right ) & =1 \end {align*}
\(c_{1}\) must be zero since \(k_{1}-\omega ^{2}m_{1}=0\) only when \(\omega =\omega _{n}\) and we assume that this is not the case here. Hence\begin {align*} c_{1} & =0\\ c_{2} & =\frac {1}{\left ( m_{1}\omega ^{2}+k_{1}\right ) } \end {align*}
Therefore \(\mathbf {x}_{p}\) becomes\[ \mathbf {x}_{p}=\begin {pmatrix} x_{1}\left ( t\right ) \\ x_{2}\left ( t\right ) \end {pmatrix} _{p}=\begin {pmatrix} \frac {1}{\left ( m_{1}\omega ^{2}+k_{1}\right ) }\sin \omega t\\ 0 \end {pmatrix} \] And the general solution becomes\begin {equation} \begin {pmatrix} x_{1}\\ x_{2}\end {pmatrix} =\begin {pmatrix} A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \\ r_{1}A_{1}^{\left ( 1\right ) }\cos \left ( \omega _{1}t+\phi _{1}\right ) +r_{2}A_{1}^{\left ( 2\right ) }\cos \left ( \omega _{2}t+\phi _{2}\right ) \end {pmatrix} +\begin {pmatrix} \frac {1}{\left ( m_{1}\omega ^{2}+k_{1}\right ) }\sin \omega t\\ 0 \end {pmatrix} \tag {8} \end {equation}
In simulink, we will directly solve the system from the original formulation\begin {align*} m_{1}\ddot {x}_{1}-k_{1}x_{2}+k_{1}x_{1} & =F\left ( t\right ) \\ m_{2}\ddot {x}_{2}+\left ( k_{1}+k_{2}\right ) x_{2}-k_{1}x_{1} & =0 \end {align*}
or\begin {align*} \ddot {x}_{1}-\frac {k_{1}}{m_{1}}x_{2}+\frac {k_{1}}{m_{1}}x_{1} & =\frac {F\left ( t\right ) }{m_{1}}\\ \ddot {x}_{2}+\frac {\left ( k_{1}+k_{2}\right ) }{m_{2}}x_{2}-\frac {k_{1}}{m_{2}}x_{1} & =0 \end {align*}
Hence\begin {align*} \ddot {x}_{1} & =\frac {F\left ( t\right ) }{m_{1}}+\frac {k_{1}}{m_{1}}x_{2}-\frac {k_{1}}{m_{1}}x_{1}\\ \ddot {x}_{2} & =-\frac {\left ( k_{1}+k_{2}\right ) }{m_{2}}x_{2}+\frac {k_{1}}{m_{2}}x_{1} \end {align*}
The simulink block diagram will be as follows for the unit step input
For an initial run with parameters \(m_{1}=m_{2}=k_{1}=k_{2}=1\) I get this warning below
EDU>> simulink Warning: Using a default value of 0.2 for maximum step size. The simulation step size will be equal to or less than this value. You can disable this diagnostic by setting 'Automatic solver parameter selection' diagnostic to 'none' in the Diagnostics page of the configuration parameters dialog.
And this is the output for \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) for the unit step response
To verify the above output from Simulink, I solved the same coupled differential equations for zero initial conditions numerically (using a numerical differential equation solver) and plotted the solution for \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) and the result matches that shown above by simulink. Here is the code the plot as a result of this verification
And the output for \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) is as follows
To verify the above output from Simulink, The same coupled differential equations were solved numerically for zero initial conditions numerically and the solution plotted for \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) and the result was found to match that shown above by simulink. Here is the code used to do the verification.
The simulink block diagram will be as follows for the \(\sin \omega t\) input
For an initial run with parameters \(m_{1}=m_{2}=k_{1}=k_{2}=1\) this is the output for \(x_{1}\left ( t\right ) \) and \(x_{2}\left ( t\right ) \) and showing the input signal at the same time
A coupled system of two masses and springs was analyzed using Simulink. The simulation was done for one set of parameters (masses and stiffness). Simulink made the simulation of this system under different loading conditions easy to do. The 2 masses response were recorded using simulink scope and the signals captured on the same plot to make it easy to compare the response of the first mass to the second mass.
The analytical analysis was more time consuming than actually making the simulation in simulink. The ability to easily change different sources to the system was useful as well as the ability to change the frequency of the input and immediately see the effect on the response.
This was my first project using Simulink, and I can see that this tool will be useful to learn more as it allows one to easily analyze engineering problems.