We want to measure the dynamical characteristics of a SDOF building system, i.e., its mass, its damping coefficient and its elastic stiffness.
To this purpose, we demonstrate that is sufficient to measure the steady-state response of the SDOF when subjected to a number of harmonic excitations with different frequencies.
The steady-state response is characterized by its amplitude, $ρ$ and the phase delay, $θ$, as in $x_{SS}(t) = ρ \sin(ωt − θ)$.
A SDOF structural system is excited by a vibrodyne that exerts a harmonic force $p(t) = p_o\sin ωt$, with $p_o = 2.224\,{}$kN at different frequencies, and we can measure the steady-state response parameters for two different input frequencies, as detailed in the following table.
$i$ | $ω_i$ (rad/s) | $ρ_i$ (μm) | $θ_i$ (deg) | $\cos θ_i$ | $\sin θ_i$ |
---|---|---|---|---|---|
1 | 16.0 | 183.0 | 15.0 | 0.9660 | 0.2590 |
2 | 25.0 | 368.0 | 55.0 | 0.5740 | 0.8190 |
We start from the equation for steady-state response amplitude,
$$\rho=\frac{p_o}{k}\frac{1}{\sqrt{(1-\beta^2)^2+(2\zeta\beta)^2}}$$where we collect $(1-\beta^2)^2$ in the radicand in the right member,
$$\rho=\frac{p_o}{k}\frac{1}{1-\beta^2}\frac{1}{\sqrt{1+[2\zeta\beta/(1-\beta^2)]^2}}$$but the equation for the phase angle, $\tan\vartheta=\frac{2\zeta\beta}{1-\beta^2}$, can be substituted in the radicand, so that, using simple trigonometric identities, we find that
$$\rho=\frac{p_o}{k}\frac{1}{1-\beta^2}\frac{1}{\sqrt{1+\tan^2\vartheta}}= \frac{p_o}{k}\frac{\cos\vartheta}{1-\beta^2}.$$With $k(1-\beta^2)=k-k\frac{\omega^2}{k/m}=k-\omega^2m$ and using a simple rearrangement, we eventually have
from scipy import matrix, sqrt, pi, cos, sin, set_printoptions
p0 = 2224.0 # converted from kN to Newton
rho1 = 183E-6 ; rho2 = 368E-6 # converted from μm to m
w1 = 16.0 ; w2 = 25.0
th1 = 15.0 ; th2 = 55.0
d2r = pi/180.
cos1 = cos(d2r*th1) ; cos2 = cos(d2r*th2)
sin1 = sin(d2r*th1) ; sin2 = sin(d2r*th2)
# the unknowns are k and m
# coefficient matrix, row i is 1, omega_i^2
coeff = matrix(((1, -w1**2),(1, -w2**2)))
# kt i.e., know term, cos(theta_i)/rho_i * p_0
kt = matrix((cos1/rho1,cos2/rho2)).T*p0
print(coeff)
print(kt)
[[ 1. -256.] [ 1. -625.]] [[ 11738901.84517425] [ 3466396.72403458]]
k_and_m = coeff.I*kt
k, m = k_and_m[0,0], k_and_m[1,0]
wn2, wn = k/m, sqrt(k/m)
print(' k m wn2 wn')
print(k, m, wn2, wn)
k m wn2 wn 17478092.3899 22418.7130654 779.620682905 27.9216883964
Using the previously established relationship for $\cos\vartheta$, we can write $\cos\vartheta=k(1-\beta^2)\frac{\rho}{p_o}$, from the equation of the phase angle (see above), we can write $\cos\vartheta = \frac{1-\beta^2}{2\zeta\beta}\sin\vartheta$, and finally
$$\frac{\rho k}{p_o}=\frac{\sin\vartheta}{2\zeta\beta}, \quad\text{hence}\quad \zeta=\frac{p_o}{\rho k}\frac{\sin\vartheta}{2\beta}$$Lets write some code that gives us our two wstimates
z1 = p0*sin1/rho1/k/2/(w1/wn)
z2 = p0*sin2/rho2/k/2/(w2/wn)
print(z1*100, z2*100)
15.7028177716 15.8171824682
We have seen that our two estimates for $\zeta$ are sligtly different, this is due to the fact that every measurement is slightly approximated...
One can partly obviate using a number of measurements larger, or even much larger, than the number of parameters s/he's trying to determine.
In our case, for each set of $M$ observations (an observation for us is a set of three values, $\omega_i, \rho_i, \theta_i$) we can write a linear equation in the $N=2$ unknowns $k$ and $m$ so that, with $b_i = p_o \cos\theta_i/\rho_i$, we can write
$$Ax-b=0$$$Ax-b=0$ is a set of $M$ equations in $N<M$ variables that cannot have an exact solution but we can introduce a vector of errors (or residuals ) with respect to an approximate solution $y \approx x$,
$$e = A\,y - b$$and assume (reasonably) that the best solution is the one that minimizes the scalar norm of the residual
$$\left|e\right|=(y^TA^T- b^T)(A\,y - b) = y^TA^TA\,y - 2 y^TA^Tb +b^Tb$$Minimizing the norm is equivalent to take the derivatives of the above expression with respect to each of the $y_i$ unknowns and equate every derivative to zero,
$$\frac{\partial |e|}{\partial y_i} = 0,\quad\text{for }i=1,\ldots,N$$and the set of equations above can be written using matrix notation:
$$A^TA\,y - A^Tb = 0.$$Solving the above equatio gives $y =(A^TA)^{-1}A^Tb$, and $y$ is the vector of parameters that minimizes the norm of the residual.
The usual name of the procedure that I have sketched is, of course, "Least Squares Parameter Estimation".
We'll deal again with the derivative of a quadratic form in the future...