Application of first-order sensitivity analysis to stabilisation of unstable continuous multi-degree-of-freedom parametric periodic systems

Wójcicki, Zbigniew

doi:10.2478/sgem-2025-0014

Full Article

1.

Introduction

The sensitivity analysis of parametric periodic systems can be an interesting theoretical problem in itself. However, the most important feature of parametric periodic systems is the instability phenomenon, which can be observed for particular values of the system parameters. Resonance vibrations in unstable parametric systems are very dangerous. Thus, stabilisation of unstable systems is usually the most important practical problem.

A full review and presentation of the existing state of knowledge in the field of parametric vibrations, with particular emphasis on stability and sensitivity analysis, will be presented in the author's next paper, entitled ‘Application of second-order sensitivity analysis to stabilisation of unstable continuous multi-degree-of-freedom parametric periodic systems’, submitted for publication in Studia Geotechnica et Mechanica together with the current work. Over a hundred papers are discussed there. However, only a few of these, which deal with the sensitivity of periodic parametric systems, are of importance from the point of view of the implementation of the goal assumed in this work.

For instance, Gu et al. [1, 2] calculate derivatives of eigenexponents, while Seyranian et al. [3] employs the sensitivity analysis of multipliers. In paper [2], the method of determination of first-order derivatives of characteristic exponents is presented. This paper contains an improvement of the method presented in [1]. This improvement allows to determine correctly the derivatives of characteristic exponents with respect to those system parameters on which the parametric excitation period depends or of which the period is itself a design parameter.

In this article, the first-order sensitivity analyses with respect to those parametric system parameters that can influence the stability/instability of the system were performed. Finally, it remains possible to determine those parameters of the system whose influence on the stabilisation procedure of such systems could be the greatest. This work continues the topics developed by the author in his earlier works, among others, in [4,5,6,7,8].

In the current paper, the original method of first-order sensitivity analysis of parametric periodic systems was formulated. There are two ways of performing stability analysis of parametric periodic systems. Both of them are connected with Floquet theory [9]. The first way is to use Lyapunov characteristic exponents, the second is to use multipliers, which are the complex eigenvalues of the monodromy matrix. The method applied in this paper is based on a sensitivity analysis of the absolute values of multipliers. From the mathematical point of view, sensitivity analysis of multipliers is the calculation of eigenderivatives with the use of derivatives of the monodromy matrix. Eigenderivatives are extremely useful for determining the sensitivities of the dynamic response to the system parameters variations. The method's innovation is the idea to achieve the sensitivity equation by analytically calculating the derivative of the homogeneous parametric equation of motion with respect to the design parameter. Then, by solving the non-homogeneous parametric sensitivity equation obtained in this way, to evaluate the first derivative of the monodromy matrix and finally the first derivatives of multipliers.

Examples of this method's implementation are also presented in this paper.

2.

Linear parametric periodic system and its stability

A linear non-homogeneous periodic parametric system of an n linear second-order differential equation of motion can be written as a first-order system (1) $\begin{matrix} \dot{x} (t) = A (t) x (t) + f (t) \\ A (t) = \begin{matrix} 0 & I \\ - B^{- 1} (t) K (t) & - B^{- 1} (t) C (t) \end{matrix}\} f (t) = \begin{matrix} 0 \\ B^{- 1} F (t) \end{matrix}\} \end{matrix}$ \matrix{{{\bf{\dot x}}(t) = {\bf{A}}(t){\bf{x}}(t) + {\bf{f}}(t)} \cr {{\bf{A}}(t) = \left\{{\matrix{{\bf{0}} & {\bf{I}} \cr {- {{\bf{B}}^{- 1}}(t){\bf{K}}(t)} & {- {{\bf{B}}^{- 1}}(t){\bf{C}}(t)} \cr}} \right\}{\bf{f}}(t) = \left\{{\matrix{{\bf{0}} \cr {{{\bf{B}}^{- 1}}{\bf{F}}(t)} \cr}} \right\}} \cr} where a 2n×2n system matrix A(t) is periodic with period T and 2n-dimensional vectors x(t) and f(t) are a vector of state and an external excitation, respectively. B(t), C(t), and K(t) are square n×n real matrices of inertia, damping, and stiffness, respectively, and F(t) is an n-dimensional excitation column vector.

Based on the Floquet theory [9], the solution of the homogeneous equation corresponding to Eq. (1)₁ has the form (2) $\begin{matrix} x (t) = X (t) x (0) & X (0) = I & X (t) = I + \int_{0}^{t} A (τ) X (τ) d τ \end{matrix}$ \matrix{{{\bf{x}}(t) = {\bf{X}}(t){\bf{x}}(0)} & {{\bf{X}}(0) = {\bf{I}}} & {{\bf{X}}(t) = {\bf{I}} + \int\limits_0^t {{\bf{A}}(\tau){\bf{X}}(\tau)d\tau}} \cr} where X(t) is a standard fundamental solution matrix or a state transition matrix normalised at zero. A steady 2n×2n monodromy matrix D is defined as the value of the fundamental matrix X(t) at the time point t = T, i.e., D = X(T).

Solving the right- and left-side eigenproblem of the monodromy matrix, i.e., (3) $\begin{matrix} \begin{matrix} (D - ρ I) r = 0 & 1^{T} (D - ρ I) = 0^{T} \end{matrix} \\ {ρ} = diag (p_{1}, p_{2}, \dots, p_{2 n}) \end{matrix}$ \matrix{{\matrix{{({\bf{D}} - \rho {\bf{I}}){\bf{r}} = {\bf{0}}} & {{{\bf{1}}^{\rm{T}}}({\bf{D}} - \rho {\bf{I}}) = {{\bf{0}}^T}} \cr}} \cr {\{{\bf{\rho}}\} = diag\left({{p_1},{p_2}, \ldots,{p_{2n}}} \right)} \cr} one can find the 2n right- and the left-side modal vectors r and l^T and the 2n multipliers, i.e., eigenvalues of the monodromy matrix D. Since the monodromy matrix D is real, non-singular, and asymmetrical, multipliers (Eq. (3)) are generally complex numbers.

The stability of the trivial solution of a homogeneous Eq. (1) depends on the absolute values of multipliers (Eq. (3)). From the point of view of practical application, a simplified system's stability/instability criterion is sufficient. If the absolute value of:

–
each multiplier is less than 1, the system is asymptotically stable in the Lyapunov sense,
–
at least one multiplier is greater than 1, the system is unstable in the Lyapunov sense.

The first and the second formula in Eq. (3) for all multipliers can be written in the form (4)

\begin{matrix} DR = R {p} & L^{T} D = {p} L^{T} \end{matrix}

\matrix{{{\bf{DR}} = {\bf{R}}\{{\bf{p}}\}} & {{{\bf{L}}^{\rm{T}}}{\bf{D}} = \{{\bf{p}}\} {{\bf{L}}^{\rm{T}}}} \cr}

where R (det R ≠ 0) is the right-side modal matrix, whose columns are right-side eigenvectors r_k of Eq. (3) corresponding to the eigenvalue ρ_k (to simplify writing, the index k is omitted henceforth). The left-side modal matrix L (det L ≠ 0) whose columns are left-side eigenvectors l_k of Eq. (3) corresponding to the eigenvalue ρ_k has to satisfy the condition (5)

L^{T} R = I i . e . L^{T} = R^{- 1}

{{\bf{L}}^{\rm{T}}}{\bf{R}} = {\bf{I}}\,\,\,{\rm{i}}.{\rm{e}}.\,\,\,{{\bf{L}}^{\rm{T}}} = {{\bf{R}}^{- 1}}

After the operation of left-side multiplication of the first formula in Eq. (4) by L^T or of right-side multiplication of the second formula in Eq. (4) by R, and considering Eq. (5), Eq. (4) can be written in a more convenient form (6)

{ρ} = L^{T} DR

\{{\bf{\rho}}\} = {{\bf{L}}^{\rm{T}}}{\bf{DR}}

3.

First-order sensitivity analysis

3.1.

First derivative of multipliers with respect to design parameter

To simplify the analysis, similarly as in the papers [1,2,3,4,5,6,7,8], only the case when the eigenvalues of the monodromy matrix are non-repeated is considered. However, in this paper, it is assumed that the multipliers can be repeated provided that Jordan's form of monodromy matrix is diagonal.

By calculating the derivative ${()}_{p}^{'} = \partial / \partial p$ \left({} \right)_p^{'} = \partial /\partial p of both sides, for example, of the first formula in Eq. (4) with respect to design parameter p, and taking into account the equality from Eq. (3), the formula for the first derivative of multipliers is received in the form (7) ${{\tilde{n}}_{p}^{'}} = L^{T} D_{p}^{'} R + L^{T} D R_{p}^{'} - L^{T} R_{p}^{'} {\tilde{n}}$ \{{\bf{\tilde n}}_p^{'}\} = {{\bf{L}}^{\rm{T}}}{\bf{D}}_p^{'}{\bf{R}} + {{\bf{L}}^{\rm{T}}}{\bf{D}}\,{\bf{R}}_p^{'} - {{\bf{L}}^{\rm{T}}}{\bf{R}}_p^{'}\{{\bf{\tilde n}}\} where the symbol ${{\tilde{n}}_{p}^{'}}$ \{{\bf{\tilde n}}_p^{'}\} denotes the diagonal matrix of the first derivatives of the monodromy matrix's multipliers with respect to the design parameter p.

3.2.

Simplifying the formula calculation of the first derivative of multipliers

From Eq. (7), it follows that to calculate the derivatives of multipliers, in addition to solving the eigenproblem Eq. (3), the derivative of the monodromy matrix $D_{p}^{'} = \partial D / \partial p$ {\bf{D}}_p^{'} = \partial {\bf{D}}/\partial p and the derivative of the right-sided eigenmatrix $R_{p}^{'} = \partial R / \partial p$ {\bf{R}}_p^{'} = \partial {\bf{R}}/\partial p must be calculated with respect to the design parameter p. In the general case, it is not possible to present analytical formulas for the calculation of both multipliers and the R and D matrices. Consequently, the formula Eq. (7), simple from the point of view of analytical operations, becomes practically useless from the numerical point of view. This formula should be treated as a matrix equation in which there are two additional unknowns $D_{p}^{'}$ {\bf{D}}_p^{'} and $R_{p}^{'}$ {\bf{R}}_p^{'} , for which we must search for numerical algorithms that allow us to achieve the goal described by Eq. (7).

It is therefore possible and necessary to further simplify Eq. (7) for the first derivative of multipliers.

For simplicity, one can transform Eq. (7), writing it as (8) $D_{p}^{'} R - R {ρ_{p}^{'}} = - D R_{p}^{'} + R_{p}^{'} {ρ}$ {\bf{D}}_p^{'}{\bf{R}} - {\bf{R}}\{{\bf{\rho}}_p^{'}\} = - {\bf{D}}\,{\bf{R}}_p^{'} + {\bf{R}}_p^{'}\{{\bf{\rho}}\} and for the k-th multiplier, as (9) $(D_{p}^{'} - ρ_{p}^{'} I) r = - (D - ρ I) r_{p}^{'}$ \left({{\bf{D}}_p^{'} - \rho _p^{'}{\bf{I}}} \right){\bf{r}} = - \left({{\bf{D}} - \rho {\bf{I}}} \right){\bf{r}}_p^{'} Multiplying Eq. (9) from the left side by the left-side eigenvector, one can obtain (10) $1^{T} (D_{p}^{'} - ρ_{p}^{'} I) r = - 1^{T} (D - ρ I) r_{p}^{'}$ {{\bf{1}}^{\rm{T}}}\left({{\bf{D}}_p^{'} - \rho _p^{'}{\bf{I}}} \right){\bf{r}} = - {{\bf{1}}^{\rm{T}}}\left({{\bf{D}} - \rho {\bf{I}}} \right){\bf{r}}_p^{'} By virtue of Eq. (3), the expression and the right side of Eq. (10) has the value 0. So, one can write this equation in the form (11) $1^{T} (D_{p}^{'} - ρ_{p}^{'} I) r = 0 or 1^{T} D_{p}^{'} r - ρ_{p}^{'} 1^{T} r = 0$ {{\bf{1}}^{\rm{T}}}\left({{\bf{D}}_p^{'} - \rho _p^{'}{\bf{I}}} \right){\bf{r}} = 0\,\,{\rm{or}}\,\,{{\bf{1}}^{\rm{T}}}{\bf{D}}_p^{'}{\bf{r}} - \rho _p^{'}{{\bf{1}}^{\rm{T}}}{\bf{r}} = 0 where by virtue of the assumption Eq. (5)₁, i.e., l^Tr = 1, one obtains (12) $ρ_{p}^{'} = 1^{T} D_{p}^{'} r or {\tilde{n}}_{p}^{'}\} = diag (L^{T} D_{p}^{'} R)$ \rho _p^{'} = {{\bf{1}}^{\rm{T}}}{\bf{D}}_p^{'}{\bf{r}}\,\,\,{\rm{or}}\,\,\,\left\{{{\bf{\tilde n}}_p^{'}} \right\} = diag\left({{{\bf{L}}^{\rm{T}}}{\bf{D}}_p^{'}{\bf{R}}} \right) instead of formula Eq. (7). Thus, the formula for calculating the derivative of multipliers becomes much more simple. Formula Eq. (12) for calculating the derivative of the multipliers of the monodromy matrix D, shows that there is no need to know two matrices $R_{p}^{'}$ {\bf{R}}_p^{'} and $D_{p}^{'}$ {\bf{D}}_p^{'} , but only one matrix $D_{p}^{'}$ {\bf{D}}_p^{'} .

3.3.

First derivative of a monodromy matrix with respect to design parameter

The first derivative of a homogeneous equation corresponding to Eq. (1) with respect to parameter p is the non-homogeneous sensitivity equation where (13) $f (t) = A_{p}^{'} (t) x (t) A_{p}^{'} = \frac{\partial A}{\partial p}$ {\bf{f}}\left(t \right) = {\bf{A}}_p^{'}\left(t \right){\bf{x}}\left(t \right)\,\,\,\,\,{\bf{A}}_p^{'} = {{\partial {\bf{A}}} \over {\partial p}} It could be proven [5,9] that the solution of this non-homogeneous sensitivity equation with zero initial conditions is the first derivative of the fundamental solution matrix with respect to the design parameter p. Ultimately, the first derivative of the fundamental and the monodromy matrix can be written as (14) $\begin{matrix} X_{p}^{'} (t) = X (t) \int_{0}^{t} X^{- 1} (τ) A_{p}^{'} (τ) X (τ) d τ \\ D_{p}^{'} = D \int_{0}^{T} X^{- 1} (τ) A_{p}^{'} (τ) X (τ) d τ \end{matrix}$ \matrix{{{\bf{X}}_p^{'}\left(t \right) = {\bf{X}}\left(t \right)\int\limits_0^t {{{\bf{X}}^{- 1}}\left(\tau \right){\bf{A}}_p^{'}\left(\tau \right){\bf{X}}\left(\tau \right)d\tau}} \cr {{\bf{D}}_p^{'} = {\bf{D}}\int\limits_0^T {{{\bf{X}}^{- 1}}\left(\tau \right){\bf{A}}_p^{'}\left(\tau \right){\bf{X}}\left(\tau \right)d\tau}} \cr} This formula can be calculated analytically or numerically, and the result may be used to calculate derivatives of multipliers in accordance with the formula Eq. (12). The same result was obtained in [3].

3.4.

The case when the period of parametric excitation depends on the design parameter

The algorithm discussed in the previous paragraph does not include a case in which the period T of the parametric excitation depends on design parameters. This is a special case. This problem was examined theoretically in the study by Gu et al. [2]. The method presented there is, unfortunately, very complicated, and the algorithm is completely useless from the point of view of the method proposed in this work. However, the results obtained in work [2] are valuable due to the possibility of their comparison with the results obtained in this paper.

The starting point to obtain the more general formula of a special case in which the period T of the parametric excitation depends on design parameters is a formula that formally describes the monodromy matrix [9] (15) $D = I + \int_{0}^{T} A (t) X (t) dt$ {\bf{D}} = {\bf{I}} + \int\limits_0^T {{\bf{A}}\left(t \right){\bf{X}}\left(t \right)dt} In the calculation of the derivative of the monodromy matrix, Eq. (15) with respect to the parameter p, it is assumed that not only matrices A = A(t,p) and X = X(t,p), but also the period T = T(p) is a function of the design parameter p. As a consequence of this assumption, the integration limit in the definite integral Eq. (15) becomes functionally dependent on the parameter p.

When calculating the derivative of the monodromy matrix, in this case, one must use the formula for the derivative of the integral with respect to the parameter [10] (16) $\frac{d}{dp} (\int_{α (p)}^{β (p)} f (t, p) dt) = \int_{α (p)}^{β (p)} \frac{\partial f (t, p)}{\partial p} dt + f (β (p), p) \frac{d β}{dp} - f (α (p), p) \frac{d α}{dp}$ {d \over {dp}}\left({\int\limits_{\alpha \left(p \right)}^{\beta \left(p \right)} {f(t,p)dt}} \right) = \int\limits_{\alpha \left(p \right)}^{\beta \left(p \right)} {{{\partial f(t,p)} \over {\partial p}}dt} + f\left({\beta (p),p} \right){{d\beta} \over {dp}} - f\left({\alpha (p),p} \right){{d\alpha} \over {dp}} In accordance with Eq. (16), one can calculate the derivative of the matrix D described by Eq. (15) (17) $\begin{array}{l} {\tilde{D}}_{p}^{'} & = \frac{d D}{dp} = \int_{0}^{T (p)} \frac{\partial}{\partial p} A (t, p) X (t, p)] dt + A (T (p), p) X (T (p), p) \frac{dT}{dp} = \\ = \int_{0}^{T (p)} [A_{p}^{'} (t, p) X (t, p) + A (t, p) X_{p}^{'} (t, p)] dt + A (T) D T_{p}^{'} \end{array}$ \matrix{{{\bf{\tilde D}}_p^{'}} \hfill & {= {{d{\bf{D}}} \over {dp}} = \int\limits_0^{T(p)} {{\partial \over {\partial p}}\left[ {{\bf{A}}\left({t,p} \right){\bf{X}}\left({t,p} \right)} \right]dt + {\bf{A}}\left({T(p),p} \right){\bf{X}}\left({T\left(p \right),p} \right){{dT} \over {dp}} =}} \hfill \cr {} \hfill & {= \int\limits_0^{T(p)} {[{\bf{A}}_p^{'}(t,p){\bf{X}}(t,p) + {\bf{A}}(t,p){\bf{X}}_p^{'}(t,p)]dt + {\bf{A}}(T){\bf{D}}T_p^{'}}} \hfill \cr} Marking (18) $D_{p}^{'} = \int_{0}^{T (p)} [A_{p}^{'} (t, p) X (t, p) + A (t, p) X_{p}^{'} (t, p)] dt$ {\bf{D}}_p^{'} = \int\limits_0^{T\left(p \right)} {[{\bf{A}}_p^{'}(t,p){\bf{X}}(t,p) + {\bf{A}}(t,p){\bf{X}}_p^{'}(t,p)]dt} and (19) ${\bar{D}}_{p}^{'} = A (T) D T_{p}^{'}$ {\bf{\bar D}}_p^{'} = {\bf{A}}(T){\bf{D}}T_p^{'} where $T_{p}^{'} = \frac{dT}{dp}$ T_p^{'} = {{dT} \over {dp}} , $A_{p}^{'} = \frac{\partial A (t, p)}{\partial p}$ {\bf{A}}_p^{'} = {{\partial {\bf{A}}(t,p)} \over {\partial p}} , $X_{p}^{'} = \frac{\partial X (t, p)}{\partial p}$ {\bf{X}}_p^{'} = {{\partial {\bf{X}}(t,p)} \over {\partial p}} , and comparing Eq. (18) and Eq. (14), one can state that Eq. (18) is visibly different from Eq. (27), however, both of them have to be equivalent. Equation (18) is the result of differentiation with respect to design parameter p of the matrix D described by Eq. (15), while Eq. (14) is a solution of the sensitivity equation achieved by differentiation of Eq. (1) with respect to the same design parameter p.

However, since there is an unknown matrix $X_{p}^{'} (t, p)$ {\bf{X}}_p^{'}(t,p) in the formula Eq. (18), this formula is practically useless. Such a method (with the use of formula Eq. (18)) of computing the derivative of a monodromy matrix can only be used if the analytical form of the fundamental matrix of solutions X(t) is known. Then, the first derivative of this matrix can be calculated analytically. Thus, in the general case, the matrix $D_{p}^{'}$ {\bf{D}}_p^{'} has to be designated using the relationship Eq. (14) instead of Eq. (18).

Finally, a more general formula than Eq. (14) for the first derivative of the monodromy matrix can be written as (20) ${\tilde{D}}_{p}^{'} = D_{p}^{'} + {\bar{D}}_{p}^{'}$ {\bf{\tilde D}}_p^{'} = {\bf{D}}_p^{'} + {\bf{\bar D}}_p^{'} It can be concluded that the formula Eq. (20) differs from the earlier formula Eq. (14) by the presence of the component ${\bar{D}}_{p}^{'}$ {\bf{\bar D}}_p^{'} , which assumes values different from zero only when the period of parametric excitation is a function of the design parameter, or is itself a design parameter. In other cases, Eq. (14) remains valid.

4.

Parametric periodic systems' stabilisation method

Based on the concept of directional derivative [10], a procedure similar to that described in [3] was used. A gradient vector has been designated for the fastest decrease of the absolute value of complex multipliers. This gradient is used to calculate the change in design parameters to make the system stable. The resulting formulas can be interpreted as an expansion of the function describing the multiplier module in Taylor series, including the first two expansion members, Eq. (31). On the other hand, it could be the first two members of the formula in work [3], where the problem was solved using the small parameter method. The algorithm was tested on the same examples that were previously analysed in works [4,5,6,7,8]. In particular, the effectiveness of the method was compared when, in addition to other system parameters, the parametric excitation period was also a design parameter. This is the fundamental difference between the algorithm presented in this work and that described in the work [3].

The possibility of a one-step exit from the area of instability was also tested.

4.1.

Gradient of the absolute value of the multiplier

Directional derivative [5] of the multiplier ρ(p₁,…,p_n) = ρ(p) as a vector function of design parameters p = [p₁,…,p_n], at the starting point specified in the parameter space with vector coordinates $p_{o} = p_{1}^{o}, \dots, p_{n}^{o}]$ {{\bf{p}}_o} = \left[ {p_1^o, \ldots,p_n^o} \right] , in the direction of any vector e = [e₁,…,e_n] can be defined as follows: (21) $ρ_{e}^{'} (p_{o}) = lim_{Δ ε \to 0} \frac{ρ (p_{o} + ε e) - ρ (p_{o})}{ε} = e^{T} ρ^{'} (p_{o}) = \sum_{k = 1}^{n} ρ_{k}^{'} e_{k}$ \rho _{\rm{e}}^{'}\left({{{\bf{p}}_o}} \right) = \mathop {\lim}\limits_{\Delta \varepsilon \to 0} {{\rho \left({{{\rm{p}}_o} + \varepsilon {\rm{e}}} \right) - \rho ({{\bf{p}}_o})} \over \varepsilon} = {{\bf{e}}^{\rm{T}}}{\bf{\rho '}}({{\bf{p}}_o}) = \sum\limits_{k = 1}^n {\rho _k^{'}{e_k}} where ɛ is a small numeric parameter, ɛe = Δp is a vector of the increments of parameters' values p = [p₁,…,p_n], a vector p^′(p_o) = grad p^′(p_o) is called a scalar function gradient ρ(p) at point p_o and indicates the direction and rate of greatest growth of this function. Thus, the vector −p^′(p) indicates, in the parameter space, the direction of the fastest decrease of the multiplier value. Vector coordinates $ρ^{'} (p_{o}) = ρ_{p_{1}}^{'}, \dots, ρ_{p_{n}}^{'}]$ {\bf{\rho '}}\left({{{\bf{p}}_o}} \right) = \left[ {\rho _{{p_1}}^{'}, \ldots,\rho _{{p_n}}^{'}} \right] are partial derivatives of the functions ρ(p) at point p_o. To simplify the index ‘_o’ in Eq. (21) was omitted (it will also be skipped in the future). The formula e^Tρ^′(p_o) means the scalar product of vectors p′(p_o) and e. Vector e^T = [e₁,…,e_n] is any vector. It is also accepted in the literature that this vector is normalised, i.e., that in the parameters space it fulfils the condition (22) $e‖ = \sqrt{e_{1}^{2} + e_{2}^{2} + \dots + e_{n}^{2}} = 1$ \left\| {\bf{e}} \right\| = \sqrt {e_1^2 + e_2^2 + \cdots + e_n^2} = 1 Due to the specificity of the problem, the vibration stabilisation algorithm will need not so much multiplier gradients as gradients of their absolute values, because their absolute values provide the stability or instability of the parametric system. A partial derivative of the absolute value of the multiplier ρ = α + iβ is defined (23) $\frac{\partial ρ|}{\partial p_{k}} = \frac{\partial \sqrt{α^{2} + β^{2}}}{\partial p_{k}} = \frac{1}{\sqrt{α^{2} + β^{2}}} (α \frac{\partial α}{\partial p_{k}} + β \frac{\partial β}{\partial p_{k}}) = \frac{Re ρ}{ρ|} g_{r_{k}}^{'} + \frac{Im ρ}{ρ|} g_{i_{k}}^{'}$ {{\partial \left| \rho \right|} \over {\partial {p_k}}} = {{\partial \sqrt {{\alpha ^2} + {\beta ^2}}} \over {\partial {p_k}}} = {1 \over {\sqrt {{\alpha ^2} + {\beta ^2}}}}\left({\alpha {{\partial \alpha} \over {\partial {p_k}}} + \beta {{\partial \beta} \over {\partial {p_k}}}} \right) = {{{\mathop{\rm Re}\nolimits} \rho} \over {\left| \rho \right|}}g_{{r_k}}^{'} + {{{\mathop{\rm Im}\nolimits} \rho} \over {\left| \rho \right|}}g_{{i_k}}^{'} where k =1,…,n. Consequently, the gradient of the absolute value of the multiplier can be written down in the form of (24) $g = grad ρ| = \frac{Re ρ}{ρ|} g_{r} + \frac{Im ρ}{ρ|} g_{u}$ {\bf{g}} = grad\left| \rho \right| = {{{\mathop{\rm Re}\nolimits} \rho} \over {\left| \rho \right|}}{{\bf{g}}_r} + {{{\mathop{\rm Im}\nolimits} \rho} \over {\left| \rho \right|}}{{\bf{g}}_u} where vector g_r is the gradient of the real part, and the vector g_u is the gradient of the imaginary part of the multiplier. According to the geometric interpretation of the scalar product of two vectors, a directional derivative is a projection of a gradient vector on the direction of the vector e. According to this interpretation, the value of the directional derivative is the greatest when its direction is the same as the direction of the gradient vector. As the possibility of reducing the value of the multiplier module is sought as soon as possible, = e = −g, which, after normalizing in accordance with Eq. (22) finally gives (25) $e = - \frac{g}{g|} = - Δ g$ {\bf{e}} = - {{\bf{g}} \over {\left| {\bf{g}} \right|}} = - \Delta {\bf{g}} Therefore, the change in the value of design parameters, leading to the fastest reduction in the absolute value of the multiplier, is determined by (26) $Δ p = ε e = - ε Δ g$ \Delta {\bf{p}} = \varepsilon {\bf{e}} = - \varepsilon \Delta {\bf{g}}

4.2.

Change of multiplier value – stabilisation procedure

One can calculate the change in the value of a multiplier using the formula (27) $Δ_{I} ρ = g_{r}^{T} Δ p + i g_{u}^{T} Δ p$ {\Delta _I}\rho = {\bf{g}}_r^T\Delta {\bf{p}} + i{\bf{g}}_u^T\Delta {\bf{p}} Using Eq. (26), one can present Eq. (27) in the form (28) $Δ_{I} ρ = Δ_{I} ρ_{r} + i Δ_{I} ρ_{u} = - ε (g_{r}^{T} Δ g + i g_{u}^{T} Δ g)$ {\Delta _I}\rho = {\Delta _I}{\rho _r} + i{\Delta _I}{\rho _u} = - \varepsilon \left({{\bf{g}}_r^{\rm{T}}\Delta {\bf{g}} + i{\bf{g}}_u^{\rm{T}}\Delta {\bf{g}}} \right) The change in the absolute value of the multiplier can be calculated based on Eq. (28) (29) $Δ_{I} ρ| = \sqrt{{(Δ_{I} ρ_{r})}^{2} + {(Δ_{I} ρ_{u})}^{2}} = ε Δ_{I} g|$ {\Delta _I}\left| \rho \right| = \sqrt {{{({\Delta _I}{\rho _r})}^2} + {{({\Delta _I}{\rho _u})}^2}} = \varepsilon {\Delta _I}\left| g \right| where the designation was adopted (30) $Δ_{I} g| = \sqrt{{(g_{r}^{T} Δ g)}^{2} + {(g_{u}^{T} Δ g)}^{2}}$ {\Delta _I}\left| g \right| = \sqrt {{{({\bf{g}}_r^{\rm{T}}\Delta {\bf{g}})}^2} + {{({\bf{g}}_u^{\rm{T}}\Delta {\bf{g}})}^2}} One can obtain an approximate multiplier value after the procedure for reducing the value of its module by adding the components described in the formula Eq. (28), i.e. (31) $ρ = ρ_{o} + Δ_{I} ρ = ρ_{r} + i ρ_{u}$ \rho = {\rho _o} + {\Delta _I}\rho = {\rho _r} + i\;{\rho _u} and the value of the module is then expressed as (32) $ρ = \sqrt{ρ_{r}^{2} + ρ_{u}^{2}}$ \rho = \sqrt {\rho _r^2 + \rho _u^2} where ρ_r = ρ_or + Δ_Iρ_r and ρ_u = ρ_ou + Δ_Iρ_u is obtained based on formula Eq. (28).

Equation (31) may, on the one hand, be interpreted as a formula corresponding to the expansion of the function describing the multiplier into the Taylor series, in which the first two members of the expansion were taken into account. On the other hand, it can be interpreted as a formula corresponding to the small parameter method, as performed in the work [3], in which the first two members of the expansion would be taken into account.

5.

Examples

5.1.

Example 1: method validation

The method presented in this paper was verified using the same example that was analysed in [2]. This example, for the parametric system, is unique. There is an analytical solution for all mathematical operations associated with the computational algorithm presented in this work. This is a great advantage of this example. It is possible to objectively verify the correctness of the theory and to determine the efficiency of the method. In addition, one can directly compare the results with those obtained in works [2, 3].

5.1.1.

Stability and sensitivity analysis

A linear parametric system described by Eq. (1) is considered, in which the system matrix (33) $A (t) = \begin{matrix} a + i (a + b) cos 2 at & - ab + ib (a + b) sin 2 at \\ \frac{a}{b} + i \frac{a + b}{b} sin 2 at & a - i (a + b) cos 2 at \end{matrix}]$ {\bf{A}}\left(t \right) = \left[ {\matrix{{a + i\left({a + b} \right)\cos 2at} & {- ab + ib\left({a + b} \right)\sin 2at} \cr {{a \over b} + i{{a + b} \over b}\sin 2at} & {a - i\left({a + b} \right)\cos 2at} \cr}} \right] is periodic with a constant period T = π/a, while a and b are system parameters. Consequently, the period T is a function of the design parameter T = T(a). The fundamental matrix of the solutions of this system, monodromy matrix, modal matrices, and first-order sensitivity analysis results have the form (34) $X (t) = \begin{matrix} cos at \cdot e^{iat} & - b sin at \cdot e^{- iat} \\ \frac{1}{b} sin at \cdot e^{iat} & cos at \cdot e^{- iat} \end{matrix}] \begin{matrix} e^{(a + ib) t} & 0 \\ 0 & e^{(a - ib) t} \end{matrix}]$ {\bf{X}}\left(t \right) = \left[ {\matrix{{\cos at \cdot {e^{iat}}} & {- b\sin at \cdot {e^{- iat}}} \cr {{1 \over b}\sin at \cdot {e^{iat}}} & {\cos at \cdot {e^{- iat}}} \cr}} \right]\left[ {\matrix{{{e^{\left({a + ib} \right)t}}} & 0 \cr 0 & {{e^{\left({a - ib} \right)t}}} \cr}} \right] (35) $D = X (T, a, b) = \begin{matrix} ρ_{1} & 0 \\ 0 & ρ_{2} \end{matrix}] = \begin{matrix} e^{(a + ib) T} & 0 \\ 0 & e^{(a - ib) T} \end{matrix}] = \begin{matrix} e^{(a + ib) \frac{π}{a}} & 0 \\ 0 & e^{(a - ib) \frac{π}{a}} \end{matrix}] = ρ\}$ {\bf{D}} = {\bf{X}}\left({T,a,b} \right) = \left[ {\matrix{{{\rho _1}} & 0 \cr 0 & {{\rho _2}} \cr}} \right] = \left[ {\matrix{{{e^{\left({a + ib} \right)T}}} & 0 \cr 0 & {{e^{\left({a - ib} \right)T}}} \cr}} \right] = \left[ {\matrix{{{e^{\left({a + ib} \right){\pi \over a}}}} & 0 \cr 0 & {{e^{\left({a - ib} \right){\pi \over a}}}} \cr}} \right] = \left\{{\bf{\rho}} \right\} (36) $R = L = I$ {\bf{R}} = {\bf{L}} = {\bf{I}} (37) $D_{a}^{'} = \begin{matrix} ρ_{1 a}^{'} & 0 \\ 0 & ρ_{2 a}^{'} \end{matrix}] = \begin{matrix} - i \frac{b}{a} T \cdot ρ_{1} & 0 \\ 0 & i \frac{b}{a} T \cdot ρ_{2} \end{matrix}]$ {\bf{D}}_a^{'} = \left[ {\matrix{{\rho _{1a}^{'}} & 0 \cr 0 & {\rho _{2a}^{'}} \cr}} \right] = \left[ {\matrix{{- i{b \over a}T \cdot {\rho _1}} & 0 \cr 0 & {i{b \over a}T \cdot {\rho _2}} \cr}} \right] (38) $D_{b}^{'} = \begin{matrix} ρ_{1 b}^{'} & 0 \\ 0 & ρ_{2 b}^{'} \end{matrix}] = \begin{matrix} - iT \cdot ρ_{1} & 0 \\ 0 & iT \cdot ρ_{2} \end{matrix}]$ {\bf{D}}_b^{'} = \left[ {\matrix{{\rho _{1b}^{'}} & 0 \cr 0 & {\rho _{2b}^{'}} \cr}} \right] = \left[ {\matrix{{- iT \cdot {\rho _1}} & 0 \cr 0 & {iT \cdot {\rho _2}} \cr}} \right] Further testing of the algorithm presented in this paper is carried out first in the case of sensitivity analysis with respect to the parameter b. This parameter does not affect the period of the parametric excitation. The results obtained in accordance with Eq. (14) should therefore be correct. The algorithm described in this formula was programmed in the Mathematica computer system. This allowed conducting analytical calculations using the symbolic procedures. The results achieved were consistent with those obtained in work [3].

The next calculations were carried out with the use of the same procedure (in accordance with the formula Eq. (14)) by studying the sensitivity with respect to the parameter a, on which the period of the parametric excitation depends. The following form of solution was obtained (39) $D_{a}^{'} = \begin{matrix} (1 + i) \frac{π}{a} e^{(1 + i \frac{b}{a}) π} & - \frac{b}{a} π e^{(1 - i \frac{b}{a}) π} \\ \frac{1}{ab} π e^{(1 + i \frac{b}{a}) π} & (1 - i) \frac{π}{a} e^{(1 - i \frac{b}{a}) π]} \end{matrix}]$ {\bf{D}}_a^{'} = \left[ {\matrix{{\left({1 + i} \right){\pi \over a}{e^{\left({1 + i{b \over a}} \right)\pi}}} & {- {b \over a}\pi {e^{\left({1 - i{b \over a}} \right)\pi}}} \cr {{1 \over {ab}}\pi {e^{\left({1 + i{b \over a}} \right)\pi}}} & {\left({1 - i} \right){\pi \over a}{e^{\left({1 - i{b \over a}} \right)\pi ]}}} \cr}} \right] This is a wrong result! The matrix $D_{a}^{'}$ {\bf{D}}_a^{'} is not even a diagonal one. One should get the result described by Eq. (37). Such an error would have to occur in the case of using the method presented in the paper [3], because the algorithm presented in this paper is not valid when the sensitivity is calculated with respect to the parameter, on which the period of the parametric excitation depends.

Ultimately, the calculations were repeated using a generalised algorithm described by the formula Eq. (20). The obtained results are identical to those achieved by analytical calculation of the derivatives, i.e. Eqs. (37) and (38).

5.1.2.

Analysis of the system stability

This example illustrates that two complex multipliers are joined by Eq. (35). Thus, their modules are the same, and a system stabilisation procedure can be carried out for any of them. It was therefore assumed that the system multiplier is (40) $ρ_{1, 2} = e^{(a \pm ib) T} = e^{aT} e^{\pm ibT} = e^{aT} (cos bT \pm i sin bT)$ {\rho _{1,2}} = {e^{\left({a \pm ib} \right)T}} = {e^{aT}}{e^{\pm ibT}} = {e^{aT}}\left({{\rm{\;cos\;}}bT \pm i{\rm{\;sin\;}}bT} \right) The absolute value of the multiplier is therefore described by the formula (41) $ρ_{1, 2}| = e^{aT}$ \left| {{\rho _{1,2}}} \right| = {e^{aT}} which functionally determines the relationship of its value to design parameters. Whether the parametric system is stable or unstable is determined, in accordance with Eq. (41), by the value of the power exponent of the number e.

If the product aT is:

–
greater than zero (aT > 0), the system will be unstable,
–
smaller than zero (aT < 0), the system will be stable.

Since the parametric excitation period in this example is T = π/a, product of aT = π > 0 and the system is unstable!

Changing the parameter b value has no effect on the stability of the solution, because the value of the multiplier module does not depend on the parameter b.

To compare the results obtained by the method of testing the values of multiplier modules with the results obtained by the study of Lyapunov exponents (as was done in the frequently cited paper [2]), one must convert one into the other.

The relationship between the multiplier ρ and characteristic exponents λ as complex numbers is described by the formula (42) $\begin{matrix} λ = \frac{1}{T} Ln ρ = \frac{1}{T} ln ρ| + i (arg ρ + 2 n π)], \\ n = 0, \pm 1, \pm 2, \dots \end{matrix}$ \matrix{{\lambda = {1 \over T}{\rm{Ln}}\;\rho = {1 \over T}\left[ {{\rm{ln}}\left| \rho \right| + i\left({{\rm{arg}}\,\rho + 2n\pi} \right)} \right],} \cr {n = 0, \pm 1, \pm 2, \ldots} \cr} Lyapunov characteristic exponents (abbreviated as Lyapunov exponent) are the real part of the characteristic exponent Eq. (42), i.e. (43) $Re λ = \frac{1}{T} \ln ρ| = \frac{1}{2 T} \ln {(Re ρ)}^{2} + {(Im ρ)}^{2}]$ {\mathop{\rm Re}\nolimits} \lambda = {1 \over T}\ln \left| \rho \right| = {1 \over {2T}}\ln \left[ {{{\left({{\mathop{\rm Re}\nolimits} \rho} \right)}^2} + {{\left({{\mathop{\rm Im}\nolimits} \rho} \right)}^2}} \right] By inserting the result of Eq. (41) into Eq. (43), the value of Lyapunov exponent is obtained, which is Re λ = a > 0. This means that the system is unstable — the same conclusion as the one reached through multipliers analysis.

A comparison of results with those obtained in the work [2] also requires the analytical calculation of the first derivatives of Lyapunov exponents (44) $λ^{'} = - \frac{1}{T^{2}} \frac{\partial T}{\partial p} \ln ρ + \frac{1}{T} \frac{1}{ρ} \frac{\partial ρ}{\partial p}$ \lambda ' = - {1 \over {{T^2}}}{{\partial T} \over {\partial p}}\ln \rho + {1 \over T}{1 \over \rho}{{\partial \rho} \over {\partial p}} Substituting the results obtained in Eq. (35) into formula (44) and taking into account that ∂T/∂a=−π/a² and ∂T/∂b = 0 finally (45) $λ_{a}^{'} = - \frac{1}{T^{2}} (- \frac{π}{a^{2}}) (a \pm ib) T + \frac{1}{T} \frac{1}{ρ} (\mp i \frac{b}{a} T ρ) = 1$ \lambda _a^{'} = - {1 \over {{T^2}}}\left({- {\pi \over {{a^2}}}} \right)\left({a \pm ib} \right)T + {1 \over T}{1 \over \rho}\left({\mp i{b \over a}T\rho} \right) = 1 (46) $λ_{b}^{'} = \frac{1}{T} \frac{1}{ρ} (\pm iT ρ) = \pm i$ \lambda _b^{'} = {1 \over T}{1 \over \rho}\left({\pm iT\rho} \right) = \pm \;i These results are consistent with those received in work [2]. At the same time, it is a confirmation that the method of testing the sensitivity of a parametric system can be implemented in two ways. The method of calculating the derivatives of the characteristic exponents of a parametric system gives the same results as the method of calculating the derivatives of multipliers.

5.1.3.

Stabilisation procedure

First, calculations were made according to the algorithm presented in the work [3]. It was assumed that there are only two design parameters in the system, i.e. p = [a,b]. The gradient determined according to Eq. (24) is then a vector with two coordinates. Since the period of the parametric excitation T = π/a is functionally dependent on the parameter a , one obtains (47) $g = {\begin{matrix} T e^{aT} & 0 \end{matrix}]}^{T} = {\begin{matrix} \frac{π e^{π}}{4} & 0 \end{matrix}]}^{T}$ {\bf{g}} = {\left[ {\matrix{{T\;{e^{aT}}} & 0 \cr}} \right]^{\rm{T}}} = {\left[ {\matrix{{{{\pi {e^\pi}} \over 4}} & 0 \cr}} \right]^{\rm{T}}} Finally, the vector Δg in Eq. (25) has a form (48) $Δ g = {[\begin{matrix} 1 & 0 \end{matrix}]}^{T}$ \Delta {\bf{g}} = {[\matrix{1 & 0 \cr} ]^{\rm{T}}} This result means that the change of the parameter a can stabilise the system. This is, of course, an incorrect result, which can easily be demonstrated. Under Eq. (41) and under the validity of the assumption that T = π/a, the absolute value of multipliers can be written as | ρ_1,2 | = e^π > 1. Therefore, the system is unstable, and the value of the multiplier module is not dependent on any design parameter. Consequently, all derived multiplier modules due to design variables are zero, and the gradient has to be a zero vector. It is therefore not possible to stabilise the system. This is contrary to the result of Eq. (48), which completes the evidence.

Once again, the cause of the error is that the algorithm presented in the work [3] does not take into account the case when the period of parametric excitation is a function of the parameter a, due to which the sensitivity of the system is being studied. In the analysed example, it is assumed that T = T(a) = π/a. The same calculations made according to the generalised algorithm presented in the work confirm the theoretical prediction that the gradient is a zero vector.

In the next calculations, it was then assumed that the vector of design variables, in addition to the earlier one, contains one more parameter, d (p = [a,b,d]). It is assumed that the period of parametric excitation depends additionally on its value, according to the formula T(a,d) = d/a, and that the initial value of the parameter d = π.

It is worth noting that performing the calculations in this variant is not possible according to the algorithm presented in the work [3]. The following detailed results were obtained (Eqs. (24) and (25)): (49) $g = {\begin{matrix} 0 & 0 & e^{d} \end{matrix}]}^{T} Δ g = {\begin{matrix} 0 & 0 & 1 \end{matrix}]}^{T}$ {\bf{g}} = {\left[ {\matrix{0 & 0 & {{e^d}} \cr}} \right]^{\rm{T}}}\;\;\;\;\;\;\;\Delta {\bf{g}} = {\left[ {\matrix{0 & 0 & 1 \cr}} \right]^{\rm{T}}} To reach the stability area in one step, the vector Δp must, by virtue of Eq. (26), fulfill the condition (50) $Δ p = {\begin{matrix} 0 & 0 & - d \end{matrix}]}^{T}$ \Delta {\bf{p}} = {\left[ {\matrix{0 & 0 & {- d} \cr}} \right]^{\rm{T}}} The final vector of design parameters takes the form of (51) $p_{1} = p_{o} + Δ p = {\begin{matrix} a & b & 0 \end{matrix}]}^{T}$ {{\bf{p}}_1} = {{\bf{p}}_o} + \Delta {\bf{p}} = {\left[ {\matrix{a & b & 0 \cr}} \right]^{\rm{T}}} With Eq. (51), it appears that in this case, the limit of the stable area will be reached only if the period of T = T(a,d) = 0/a = 0. Mathematically, this is the correct result, but from the point of view of physical interpretation, this means the period of parametric excitation disappears. Stabilisation of the system is therefore possible but requires complete elimination of parametric excitation. This means the de facto transition to a system with constant coefficients.

5.2.

Example 2: stabilisation of an unstable parametric system by a dynamic absorber

The problem of optimal tuning of the dynamic eliminator in a single-degree-of-freedom system with constant coefficients was formulated and solved in [11]. It was shown that if harmonic excitation is the only factor that excites the constant parameter system, the additional mass (and thus the added second degree-of-freedom) attached to the primary mass can effectively act as a dynamic eliminator of resonant vibrations of the primary mass.

Following this phenomenon, an attempt was made to stabilise the unstable parametrically excited system using an additional active (classical) dynamic absorber. ‘Active’ denotes an absorber whose parameters can be changed during the stabilisation process.

It was initially assumed that the design parameters of the system were the characteristics of the absorber only, i.e. the stiffness and damping of the absorber. It was considered that in the case of engineering structures, it is much more difficult to change the features of the system than the parameters of the vibration absorber. The adverse effects of parametric resonant vibrations often manifest themselves only after the construction of the structure. Therefore, adding an absorber is a more viable option than changing an existing structure.

In the case of the active parametric dynamic absorber, it was assumed that the frequency of parametric excitation of the system and the frequency of the absorber characteristics are the same. This time, ‘active’ denotes an absorber whose parameters can be changed during the stabilisation and/or operation process.

The effect of the length of the parametric excitation period change, which is then both a parameter of the system and the parametric absorber, was additionally studied.

This example aims to seek an answer to the question: Is a dynamic absorber able to effectively stabilise an unstable parametric system?

The tool used to achieve the goal is the method presented in this paper. First, the sensitivity of the eigenproblem of the monodromy matrix of the parametric system is analysed, and the derivatives of the modulus of the largest multiplier “responsible’ for the instability of the system are determined. Then, using the presented gradient method, the absorber parameters are automatically modified to reduce the absolute value of the largest multiplier in the fastest way. If the modulus values of all the multipliers thus become less than one, the system will leave the region of instability and become stable.

If it is not possible to stabilise the system with the use of a dynamic vibration absorber, the task can be extended to study the influence of the system parameters as well.

The effectiveness of two types of dynamic vibration absorbers is analysed: classic and parametric ones.

Symbolic and numerical calculations are performed using a computer system, Mathematica, Wolfram [12].

5.2.1

Model of the system-absorber

By virtue of the theory presented in the paper, a body with a concentrated mass M(t) (primary mass of the parametric system) with a dynamic vibration absorber (mass m) attached to it, is analysed. Thus, a two-mass damped parametric system with two dynamic degrees of freedom, q₁ and q₂, is created. The system is shown in Fig. 1.

The additional, upper mass m in Fig. 1 is the mass of the dynamic vibration absorber connected in parallel by an elastic spring k(t) and damping bond c(t) with the primary mass M(t) (lower mass) of the parametric system, attached to the foundation by an elastic spring K(t) and damping bond C(t). The system is parametrically excited—elastic spring stiffness K(t) and damping characteristic C(t) periodically change in time with a parametric excitation frequency n_o. It is also assumed that, in a general case, the stiffness k(t) and damping characteristics c(t), connecting the mass of the vibration eliminator m with the primary mass M(t), can harmonically vary in time with the frequency of parametric excitation.

The matrix equation of motion of a parametric homogenous system shown in Fig. 1 can be generally written in the form (52) $B (t) \ddot{q} + C (t) \dot{q} + K (t) q = 0$ {\bf{B}}\left(t \right){\bf{\ddot q}} + {\bf{C}}\left(t \right){\bf{\dot q}} + {\bf{K}}\left(t \right){\bf{q}} = {\bf{0}} where B(t), C(t) and K(t) are inertia, damping and stiffness matrices and $q = {[\begin{matrix} q_{1} & q_{2} \end{matrix}]}^{T}$ {\bf{q}} = {[\matrix{{{q_1}} & {{q_2}} \cr} ]^T} is a general coordinates vector, Fig. 1.

In this example, without losing the generality of Eq. (52), it was assumed that the system characteristics of elements in Fig. 1 are described by the following values [4]:

Primary parametric system (lower mass M)

the basic mass	M = 3 × 10⁴ kg
the constant part of the stiffness of the elastic bond	K_o = 4 × 10⁸ N/m
the characteristic of the damping bond	C = 4.4 × 10⁴ Ns/m
Parametric excitation frequency Eliminator (upper mass m)	v₀−233 rad/s
the eliminator mass (1/20 of the mass M)	m = 1.5 × 10³ kg
the spring stiffens	k = 1.77 × 10⁷ N/m
the characteristic of the damping bond	c_o = 4.4 × 10⁴ N/m

5.2.2.

The classic vibration absorber

The matrix coefficients of the equation of motion Eq. (52) can be then described by the following formulas and matrix coefficients: (53) $B (t) = \begin{matrix} m & 0 \\ 0 & M \end{matrix}]$ {\bf{B}}\left(t \right) = \left[ {\matrix{m & 0 \cr 0 & M \cr}} \right] (54) $C (t) = \begin{matrix} c_{o} & - c_{o} \\ - c_{o} & c_{o} + C \end{matrix}]$ {\bf{C}}\left(t \right) = \left[ {\matrix{{{c_o}} & {- {c_o}} \cr {- {c_o}} & {{c_o} + C} \cr}} \right] (55) $K (t) = \begin{matrix} k & - k \\ - k & k + K_{o} \end{matrix}] + \begin{matrix} 0 & 0 \\ 0 & K_{1} \end{matrix}] cos (ν_{o} t)$ {\bf{K}}\left(t \right) = \left[ {\matrix{k & {- k} \cr {- k} & {k + {K_o}} \cr}} \right] + \left[ {\matrix{0 & 0 \cr 0 & {{K_1}} \cr}} \right]{\rm{\;}}\;\;{\rm{cos\;}}\left({{\nu _o}\;t} \right) where (56) $K (t) = K_{o} + K_{1} cos (ν_{o} t), (K_{o} > K_{1} = 0.5 K_{o})$ K\left(t \right) = {K_o} + {K_1}{\rm{\;cos\;}}\left({{\nu _o}\;t} \right),\;\;\;\;\;({K_o} > {K_1} = 0.5{K_o}) In order to use the theory presented in paper, the equation of motion Eq. (52) is reduced to a system of the first-order homogenous differential equations corresponding to Eq. (1)₁.

According to the theory presented in paper, the system matrix A(t) in Eq. (1)₂ has than the form (57) $A (t) = \begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ - \frac{k}{m} & \frac{k}{m} & - \frac{c_{o}}{m} & \frac{c_{o}}{m} \\ \frac{k}{M} & - \frac{(K_{o} + K_{1} cos v_{o} t) + k}{M} & \frac{c_{o}}{M} & - \frac{C + c_{o}}{M} \end{matrix}]$ {\bf{A}}\left(t \right) = \left[ {\matrix{0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr {- {k \over m}} & {{k \over m}} & {- {{{c_o}} \over m}} & {{{{c_o}} \over m}} \cr {{k \over M}} & {- {{\left({{K_o} + {K_1}{\rm{\;cos\;}}{v_o}t} \right) + k} \over M}} & {{{{c_o}} \over M}} & {- {{C + {c_o}} \over M}} \cr}} \right] On the basis of stability analysis, it was found that the initial system is unstable. This is evident from the values of the multipliers on the major diagonal of the Jordanian form of the monodromy matrix. The first multiplier is modularly larger than one (58) $J = \begin{matrix} - 1.1938 & 0 & 0 & 0 \\ 0 & - 0.66404 + 0.22579 i & 0 & 0 \\ 0 & 0 & - 0.66404 - 0.25579 i & 0 \\ 0 & 0 & 0 & - 0.69295 \end{matrix}]$ {\bf{J}} = \left[ {\matrix{{- 1.1938} & 0 & 0 & 0 \cr 0 & {- 0.66404 + 0.22579\;i} & 0 & 0 \cr 0 & 0 & {- 0.66404 - 0.25579\;i} & 0 \cr 0 & 0 & 0 & {- 0.69295} \cr}} \right] To stabilise the system, the method of first-order sensitivity analysis presented in this paper was performed. The design parameter vector is assumed to be in form (59) $p = {\begin{matrix} k & c_{o} \end{matrix}]}^{T}$ The values of the design parameter vector in start point are the values (60) $p_{0} = {\begin{matrix} 1.77 \cdot 10^{7} & 44000 \end{matrix}]}^{T}$ {{\bf{p}}_0} = {\left[ {\matrix{{1.77 \cdot {{10}^7}} & {44000} \cr}} \right]^T} The following normalised gradient vector Eq. (25) were obtained at the beginning (61) $Δ g = {\begin{matrix} - 0.0091117 & 0.99996 \end{matrix}]}^{T}$ \Delta {\bf{g}} = {\left[ {\matrix{{- 0.0091117} & {0.99996} \cr}} \right]^T} and the end of the stabilisation process (62) $Δ g = {\begin{matrix} - 0.019800 & 0.99980 \end{matrix}]}^{T}$ \Delta {\bf{g}} = {\left[ {\matrix{{- 0.019800} & {0.99980} \cr}} \right]^T} As one can see, the change in the stiffness of the eliminator has a much smaller (practically insignificant) effect on the stabilisation of the system, despite the fact that the absolute value of the element of the gradient vector corresponding to the stiffness increased during the stabilisation process. The system was able to stabilise at the level of (63) $p = {\begin{matrix} 1.77 \cdot 10^{7} & 22213 \end{matrix}]}^{T}$ {\bf{p}} = {\left[ {\matrix{{1.77 \cdot {{10}^7}} & {22213} \cr}} \right]^T} The stiffness value has hardly changed, while the eliminator damping has decreased almost twice compared to the initial value. Such a tuned eliminator guarantees the stability of the entire system, as evidenced by the values of the multipliers (64) $J = \begin{matrix} - 0.9917 - 0.02384 i & 0 & 0 & 0 \\ 0 & - 0.9917 - 0.02384 i & 0 & 0 \\ 0 & 0 & - 0.7598 + 0.2601 i & 0 \\ 0 & 0 & 0 & - 0.7598 - 0.2601 i \end{matrix}]$ {\bf{J}} = \left[ {\matrix{{- 0.9917 - 0.02384\;i} & 0 & 0 & 0 \cr 0 & {- 0.9917 - 0.02384\;i} & 0 & 0 \cr 0 & 0 & {- 0.7598 + 0.2601\;i} & 0 \cr 0 & 0 & 0 & {- 0.7598 - 0.2601\;i} \cr}} \right] The largest modulus of the multiplier has the value 0.992, which is less than 1. This indicates the stability of the system.

5.2.3.

The parametric vibration absorber

The same unstable parametric system, which was stabilised by a classical vibration absorber, was subjected to the stabilisation procedure described in this work, but this time, the characteristic of the damping bond of the vibration absorber was modified. It was assumed that (65) $c (t) = c_{o} + c_{1} cos (ν_{o} t), (c_{o} > c_{1} = 0.5 c_{o})$ c\left(t \right) = {c_o} + {c_1}{\rm{\;cos\;}}\left({{\nu _o}t} \right),\;\;\;\;\;({c_o} > {c_1} = 0.5{c_o}) Because the characteristic of the bond periodically changes in time with a frequency that equals the parametric excitation frequency, the absorber has become a parametric vibration absorber. Damping matrix in equation of motion Eq. (52) can be written as (66) $C (t) = \begin{matrix} c_{o} & - c_{o} \\ - c_{o} & c_{o} + C \end{matrix}] + \begin{matrix} c_{1} & - c_{1} \\ - c_{1} & c_{1} \end{matrix}] cos (ν_{o} t)$ {\bf{C}}\left(t \right) = \left[ {\matrix{{{c_o}} & {- {c_o}} \cr {- {c_o}} & {{c_o} + C} \cr}} \right] + \left[ {\matrix{{{c_1}} & {- {c_1}} \cr {- {c_1}} & {{c_1}} \cr}} \right]{\rm{\;cos\;}}\left({{\nu _o}\;t} \right) System matrix A(t) Eq. (57) must be modified, and now has the form (67) $A (t) = \begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ - \frac{k}{m} & \frac{k}{m} & - \frac{c_{o} + c_{1} cos ν_{o} t}{m} & \frac{c_{o} + c_{1} cos ν_{o} t}{m} \\ \frac{k}{M} & - \frac{(K_{o} + K_{1} cos ν_{o} t) + k}{M} & \frac{c_{o} + c_{1} cos ν_{o} t}{M} & - \frac{(c_{o} + c_{1} cos ν_{o} t) + C}{M} \end{matrix}]$ {\bf{A}}\left(t \right) = \left[ {\matrix{0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 1 \cr {- {k \over m}} & {{k \over m}} & {- {{{c_o} + {c_1}{\rm{\;cos\;}}{\nu _o}t} \over m}} & {{{{c_o} + {c_1}{\rm{\;cos\;}}{\nu _o}t} \over m}} \cr {{k \over M}} & {- {{\left({{K_o} + {K_1}{\rm{\;cos\;}}{\nu _o}t} \right) + k} \over M}} & {{{{c_o} + {c_1}{\rm{\;cos\;}}{\nu _o}t} \over M}} & {- {{\left({{c_o} + {c_1}{\rm{\;cos\;}}{\nu _o}t} \right) + C} \over M}} \cr}} \right] The initial values of the multipliers have hardly changed, compared to the previously described relation, Eq. (58). (68) $J = \begin{matrix} - 1.1912 & 0 & 0 & 0 \\ 0 & - 0.66352 + 0.24280 i & 0 & 0 \\ 0 & 0 & - 0.66351 - 0.24280 i & 0 \\ 0 & 0 & 0 & - 0.70447 \end{matrix}]$ {\bf{J}} = \left[ {\matrix{{- 1.1912} & 0 & 0 & 0 \cr 0 & {- 0.66352 + 0.24280\;i} & 0 & 0 \cr 0 & 0 & {- 0.66351 - 0.24280\;i} & 0 \cr 0 & 0 & 0 & {- 0.70447} \cr}} \right] Similar to the coordinates values of a normalised gradient vector, Eq. (61) (69) $Δ g = {\begin{matrix} - 0.0090000 & 0.99996 \end{matrix}]}^{T}$ \Delta {\bf{g}} = {\left[ {\matrix{{- 0.0090000} & {0.99996} \cr}} \right]^T} However, it was not possible to stabilise the system. As in the case described above, the coordinate values of the normalised gradient did not change much in the stabilisation process. However, the gradient before the normalisation had practically zero values, which indicates that the local extreme had been reached. The largest modulus has the value 1.00246 and is greater than one, which means that the system is still unstable. The limit of the area of stability was practically reached, but it was not exceeded, and, as mentioned earlier, the system remained unstable.

The automatic stabilisation test was repeated assuming a decrease in the value of the parametric part of the eliminator damping bond characteristic, i.e. it is assumed that c₁ = 0.3c_o instead of c₁ = 0.5c_o. This time, the stabilisation process was successful and, as in the case of the classic vibration eliminator, a stable system was obtained, in which the maximum value of the multiplier modulus was 0.99737. The vector of design parameters also did not differ much from that obtained for a classic vibration eliminator, i.e. (70) $p = {\begin{matrix} 1.77002 \cdot 10^{7} & 23358 \end{matrix}]}^{T}$ {\bf{p}} = {\left[ {\matrix{{1.77002 \cdot {{10}^7}} & {23358} \cr}} \right]^T} At the end, one more variant of calculations was carried out. It was assumed that all the parameters of the vibration absorber and the parametric system are design parameters, i.e. that the vector of the design parameters consisted of seven coordinates (71) $p = {\begin{matrix} k & K_{o} & K_{1} & c_{o} & C & c_{1} & T \end{matrix}]}^{T}$ {\bf{p}} = {\left[ {\matrix{k & {{K_o}} & {{K_1}} & {{c_o}} & C & {{c_1}} & T \cr}} \right]^{\rm{T}}} whose values are (72) $p_{o} = {\begin{matrix} 1.77 \cdot 10^{7} & 4 \cdot 10^{8} & 2 \cdot 10^{8} & 44000 & 44000 & 22000 & 0.026995 \end{matrix}]}^{T}$ {{\bf{p}}_o} = {\left[ {\matrix{{1.77 \cdot {{10}^7}} & {4 \cdot {{10}^8}} & {2 \cdot {{10}^8}} & {44000} & {44000} & {22000} & {0.026995} \cr}} \right]^{\rm{T}}} The last parameter is the parametric excitation period of the system T = 2π/v_o = 0.0269665 s. To perform this variant of calculations, it was necessary to use generalised formulas for calculating the derivatives of the monodromy matrix Eq. (19).

The design parameters vector after stabilizing the system is shown below: (73) $p = {\begin{matrix} 1.77 \cdot 10^{7} & 4 \cdot 10^{8} & 2 \cdot 10^{8} & 44000 & 44000 & 22000 & 0.024144 \end{matrix}]}^{T}$ {\bf{p}} = {\left[ {\matrix{{1.77 \cdot {{10}^7}} & {4 \cdot {{10}^8}} & {2 \cdot {{10}^8}} & {44000} & {44000} & {22000} & {0.024144} \cr}} \right]^{\rm{T}}} Values of the vector p differ from the values of the vector p_o practically only at the last position—corresponding to the parameter, which is the parametric excitation period. Compared to the initial value, the parametric excitation period decreased by only about 10%. Thus, in this case, a small change in a sole parameter guarantees the stability of the system. The stabilisation process was carried out instantly, which is important if the stabilisation procedure is to be carried out in practice, in which case it must take place in real time.

A normalised gradient vector contains virtually only one non-zero element (74) $Δ g = {\begin{matrix} - 3.451 \cdot 10^{- 11} & 2.451 \cdot 10^{- 11} & 1.312 \cdot 10^{- 11} & - 2.877 \cdot 10^{- 8} & - 1.725 \cdot 10^{- 9} & 7.448 \cdot 10^{- 9} & 1.000 \end{matrix}]}^{T}$ \Delta {\bf{g}} = {\left[ {\matrix{{- 3.451 \cdot {{10}^{- 11}}} & {2.451 \cdot {{10}^{- 11}}} & {1.312 \cdot {{10}^{- 11}}} & {- 2.877 \cdot {{10}^{- 8}}} & {- 1.725 \cdot {{10}^{- 9}}} & {7.448 \cdot {{10}^{- 9}}} & {1.000} \cr}} \right]^{\rm{T}}} This confirms the earlier conclusion that the only important parameter on which the rate of stabilisation practically depends is the parametric excitation period. The influence of other parameters is smaller by at least eight orders.

The Jordan form of the monodromy matrix after stabilisation of the system supports the conclusion of its stability. This is because one gets (75) $J = \begin{matrix} - 0.0017160 & 0 & 0 & 0 \\ 0 & - 0.65648 + 0.49881 i & 0 & 0 \\ 0 & 0 & - 0.65648 - 0.49881 i & 0 \\ 0 & 0 & 0 & - 0.70099 \end{matrix}]$ {\bf{J}} = \left[ {\matrix{{- 0.0017160} & 0 & 0 & 0 \cr 0 & {- 0.65648 + 0.49881\;i} & 0 & 0 \cr 0 & 0 & {- 0.65648 - 0.49881\;i} & 0 \cr 0 & 0 & 0 & {- 0.70099} \cr}} \right] The largest absolute value of the multiplier is 0.99172 and is less than 1, which indicates the stability of the system, while the absolute value of the largest multiplier could be reduced even further.

6.

Conclusions

The aim of the study was to propose an original method for stabilizing an unstable multi-degree-of-freedom parametric system. The method's innovation is the idea to achieve the non-homogeneous parametric sensitivity equation by evaluating analytically the derivative of the homogeneous parametric differential equation of motion with respect to the design parameter. Then, by solving the sensitivity equation obtained in this way, to evaluate the first derivative of the monodromy matrix and finally the first derivatives of multipliers. Ultimately, this method is based on the sensitivity analysis of absolute values of multipliers. As a result, numerical rather than analytical procedures can be used, which is a significant improvement in the case of parametric systems, in which analytical solutions practically do not exist. This procedure is based on the concept of first-order sensitivity analysis, which allows to determine those design parameters that have the greatest influence on the response of the parametric system. Next, using the gradient method, the values of selected design parameters are changed in such a way as to stabilise the parametric system as quickly as possible.

The method was modified in such a way that, in particular, it becomes possible to use the parametric excitation period also as a design parameter. The method presented in work [3] does not provide such a possibility, and in the case of dependence of the parametric excitation period on the design parameter, an error is generated.

Two basic assumptions were made in the paper:

–
the parametric systems are limited to the case of a linear parametric system in which the variability in time of its parameters is described by a continuous periodic function of time,
–
the first-order sensitivity analysis is applied only.

Two examples were analysed in the paper. This first example – method validation (the same example that was analysed in [2]), is absolutely unique for the parametric systems, since an analytical solution for all mathematical operations associated with the presented algorithm exists. This is a great advantage of this example. It is possible to objectively verify the correctness of the theory and to determine the efficiency of the method. In addition, one can directly compare the results with those obtained at works [2,3]. The comparison was realised in parallel in three ways:

–
analytical calculations performed by a human,
–
analytical calculations performed using symbolic procedures of the Mathematica system,
–
calculations performed using numerical procedures of the Mathematica system.

The consistency of the results obtained by these three approaches irrefutably proves the correctness of the method described in the paper.

The goal of the second example – stabilisation of an unstable parametric system by a dynamic absorber – was different. It is an attempt to show the possibility of the practical application of the proposed method.

In this example, the effectiveness of two types of absorbers was analysed: classical and parametric ones.

It was found that both classical and parametric vibration eliminators could be an effective tool in achieving this goal. The research shows the following additional detailed conclusions:

–
a classic eliminator can be more effective than a parametric one in the process of stabilizing parametric vibrations,
–
in the stabilisation process of a parametric system, the parametric excitation period is the parameter whose change has the greatest influence on the speed and efficiency of the stabilisation process; unfortunately, the period of parametric variability of the vibration absorber characteristic is not only a parameter of the absorber but also a parameter of the system.

The authors' next paper, submitted for publication together with the current one, consists of an extension of the method presented in this paper with second-order sensitivity analysis algorithms. Further works on the sensitivity analysis of discontinuous multi-degree-of-freedom parametric periodic systems are now in the final phase of preparations.

Application of first-order sensitivity analysis to stabilisation of unstable continuous multi-degree-of-freedom parametric periodic systems

Full Article

Paradigm

My account