# 3-dimensional flutter kinematic structural stability

Jean Lerbet, Gaëtan Hello, Noël Challamel, François Nicot, Félix Darve

## To cite this version:

Jean Lerbet, Gaëtan Hello, Noël Challamel, François Nicot, Félix Darve. 3-dimensional flutter kinematic structural stability. Nonlinear Analysis: Real World Applications, 2016, 29, pp.19-37. 10.1016/j.nonrwa.2015.10.006. hal-01321946

HAL Id: hal-01321946

https://hal.science/hal-01321946

Submitted on 6 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3-dimensional flutter kinematic structural stability

J. Lerbet$^{a,*}$, G. Hello$^b$, N. Challamel$^c$, F. Nicot$^d$, F. Darve$^e$

$^a$ IBISC, UFRST-UEVE, 40, rue du Pelvoux CE 1455, 91020 Evry Courcouronnes cedex, France

$^b$ LMEE, UFRST-UEVE, 40, rue du Pelvoux CE 1455, 91020 Evry Courcouronnes cedex, France

$^c$ Université Européenne de Bretagne Université de Bretagne Sud LIMATB-UBS-Lorient Centre de Recherche Rue de Saint Maudé-BP 92116, 56321 Lorient cedex, France

$^d$ IRSTEA, ETNA-Geomechanics Group, 2, rue de la papeterie, 38042 St Martin d’Heres cedex, France

$^e$ Univ. Grenoble Alpes, 3SR, F-38000 Grenoble, France

Having recalled the kinematic structural stability (ki.s.s) issue and its solution for divergence-type instability, we address the same problem for flutter-type instability for the minimal required configuration of dimensions—meaning 3 degree of freedom systems. We first get a sufficient non optimal condition. In a second time, the complete issue is tackled by two different ways leading to same results. A first way using calculations on Grassmann and Stiefel manifolds that may be generalized for any dimensional configuration. A second way using the specific dimensional configuration is brought back to calculations on the sphere. Differences with divergence ki.s.s are highlighted and examples illustrate the results.

## 0. Introduction

This paper deals with the so-called kinematic structural stability (ki.s.s.) for the flutter of non conservative elastic discrete systems. In a previous recent paper (see [1]), the ki.s.s. problematic was formulated in its generality and the solution for the divergence criterion for conservative as well as for non conservative elastic discrete systems has also been given by use of two independent ways. The first one has been proposed for some years by using the formula of Schur’s complements, using Lagrange multipliers for introducing the kinematic constraints (see for example [2,3]). The second approach [1] is based upon a variational formulation of the divergence criterion and the explicit elimination of Lagrange’s multipliers associated to the additional kinematic constraints. Both approaches lead (fortunately!) to the same results: for conservative elastic systems, the ki.s.s. is universal (as it was for long time known) and can be proved by the use of Rayleigh’s quotient and Courant’s Minimax results: in fact, adding a kinematic constraint on a conservative

<footnote>* Corresponding author. Tel.: +33 1 69 47 75 69; fax: +33 1 69 47 75 03.
E-mail address: jean.lerbet@ibisc.univ-evry.fr (J. Lerbet).</footnote>
system cannot destabilize an equilibrium position. We may first remark that, in this last phrase as in the
whole paper, possible gyroscopic effects (and especially usual stabilization effects of gyroscopic forces) are
not taken into account (see [4]). If not, it is then easy to exhibit a counter-example to the sentence about this
ki.s.s. property. We may also mention that Tarnai clearly showed (see [5] for example) that this property fails
even for conservative systems if additional kinematic constraints change the considered equilibrium position.
Kinematic structural stability refers, by definition, to a given equilibrium position of a mechanical system
Σ that must not be changed by adding kinematic constraints: only the eventual change of stability of this
equilibrium configuration is investigated. On the contrary, as already mentioned by Thompson [6] in 1982
but never systematically investigated before a set of recent papers ([2,3,1] for example), the non universal
divergence-type ki.s.s. is characteristic of the nonconservativity of Σ and the main result reads so:
for non conservative systems, the divergence-type ki.s.s. (or more concisely the divergence ki.s.s.) is only
conditional according to the second order work criterion: as long as the symmetric part of the stiffness matrix
remains definite positive, the ki.s.s. holds and no additional kinematic constraint may destabilize the system
by divergence. As soon as the isotropic cone is not nil, the invert image by the stiffness matrix of a vector
chosen on this cone provides a constraint that destabilizes the system by divergence.

In this paper, we focus on the flutter criterion. Because this issue is much more difficult, we only deal
meanwhile with 3 dof systems subjected to one additional kinematic constraint and because flutter may occur
only for at least 2 dof systems, this configuration of dimensions is the minimal one required to question the
ki.s.s issue for the flutter instability. Contrary to the divergence ki.s.s issue, we did not find out a pure
algebra-geometric reasoning allowing to solve the problem and differential calculations are definitively
necessary. The flutter ki.s.s. is brought back to an optimization problem with vector subspaces as
optimization variables. That leads us to use differential geometry tools as Grassmann manifolds even if
a parametrization through the sphere of the 3 dimensional euclidean space may be used remaining careful
that a same two dimensional vector space has two unit normal vectors. Obviously this parametrization could
not be used in higher dimensions whereas the reasoning with Grassmann and Stiefel manifolds is more easily
generalizable. For precisions on Grassmann and Stiefel manifolds especially for applications to numerical
methods and optimization issues see for example [7,8].

Calculations are done here by both methods and show that the flutter ki.s.s. is neither universal nor con-
ditional but must be handled case by case. There are systems Σ = Σ<sub>free</sub> where all the associated constrained
systems Σ<sub>C</sub> are more stable than the initial free system meaning that the critical flutter load value p<sup>*</sup><sub>fl</sub> for
the free system Σ<sub>free</sub> is lower than the critical flutter load value p<sup>*</sup><sub>fl,C</sub> for any constrained system Σ<sub>C</sub>. On the
contrary, there are systems Σ where at least one associated constrained system Σ<sub>C</sub> is less stable than the
initial free system meaning that the critical flutter load value p<sup>*</sup><sub>fl</sub> for the free system Σ is higher than the
critical flutter load value p<sup>*</sup><sub>fl,C</sub> for the considered constrained system Σ<sub>C</sub>.

The paper is organized as follows. In Section 1, the general ki.s.s. problematic and its solution for
divergence type instabilities are first quickly summarized. Then the flutter ki.s.s. issue is formalized leading
to an optimization problem for a well-defined function Φ on a Grassmann manifold. In Section 2, some
calculations, used subsequently, are done which lead to a significant sufficient condition for preserving
flutter ki.s.s. This algebraic condition involves spectrum of both symmetric and skew symmetric parts of the
operator. In Section 3, the main calculations are done leading the critical points of Φ. As mentioned above,
two ways are used that lead to the same results. The first one is more condensed but uses less known tools
of differential geometry. It may be generalized to higher dimensional issues. The second way is more usual
by use of differential calculations on the sphere but it leads to more complicated calculations and cannot be
generalized for flutter ki.s.s. issue in higher dimension. The general results show again that the flutter ki.s.s.
is controlled through a competition between the symmetric and skew symmetric parts of a single operator
expressed with respect to the stiffness and the mass matrices of the system. In the fourth and last section,
numerical calculations for the paradigmatic 3 dof Ziegler system illustrate the general analytic results.
# 1. K.I.S.S.

## 1.1. Divergence ki.s.s

Let $\Sigma = \Sigma_{free}$ be a mechanical system, $q = (q_1, ..., q_n)$ a (local) coordinate system of its configuration space, $p$ a loading parameter and $q^e$ an equilibrium configuration of $\Sigma$ subjected to a load system characterized to simplify by a single dimensionless parameter $p$. In this paper, only linear stability is investigated. Suppose then that $q^e$ is linearly stable in the Lyapunov sense. That means that the dynamical system

$$M \ddot{X} + K(p)X = 0 \quad (1)$$

is Lyapunov stable where $M$ is the symmetric positive definite mass matrix and $K(p)$ the stiffness matrix is without any property and especially generically non symmetric. By using the square root $S$ of $M$ we convert (1) into

$$\ddot{Y} + \tilde{K}(p)Y = 0 \quad (2)$$

where $\tilde{K}(p) = S^{-1}K(p)S^{-1}$ contains in a unique matrix the whole dynamics. The load parameter $p$ is suppose monotone increasing and often we suppose that for $p = 0$ the system is conservative stable implying that $K(0)$ is symmetric positive definite. If the system becomes unstable for $p = p^*$, the domain of stability is then $[0, p^*[$ with eventually $p^* = +\infty$.

For divergence stability the equation becomes

$$K(p)X = 0 \quad \text{or equivalently} \quad \tilde{K}(p)Y = 0. \quad (3)$$

A family of $\ell$ linear kinematic constraints is a family $C = (C_1, ..., C_\ell)$ of linear forms identified by the scalar product with column vectors so that these constraints are equivalent to a matrix $C = \text{mat}(C_1...C_\ell) \in \mathcal{M}_{n\ell}(\mathbb{R})$. The constrained system will be denoted by $\Sigma_C$ and as mentioned in the introduction, $0$ is supposed to be still an equilibrium position of $\Sigma_C$. Analogously, $[0, p_C^*[$ is the stability interval of the configuration $0$ for $\Sigma_C$. The kinematic structural stability (ki.s.s.) refers to the preserving of the stability of any system $\Sigma_C$ produced from $\Sigma$ by adding any family of kinematic constraints $C$ to $\Sigma$. That means that for the same value $p$ of the loading, $q_e$ is still an equilibrium configuration of the constrained system $\Sigma_C$ and that $q^e$ is still Lyapunov stable as equilibrium configuration of $\Sigma_C$.

For elastic conservative systems, Courant Minimax results about Rayleigh's quotient may then be easily translated as universal ki.s.s. (for divergence and for flutter as well because the instability may only occur by divergence!): $p^* \le p_C^* \forall C$.

On the contrary, for non conservative elastic systems, the ki.s.s issue is more complicated because, at least for linear stability, two modes of instability may occur: divergence instability and flutter instability. We then separate the ki.s.s into two types of ki.s.s.: one for divergence and one other for flutter.

The recent results mentioned in the introduction show that the divergence ki.s.s. is conditional, the condition being that the symmetric part of the stiffness matrix remains positive definite: it is nothing else that the second order work criterion that then appears as the optimal criterion for divergence ki.s.s. (see [2,3,1]). Denoting by $p_{div}^*$ the critical divergence load of $\Sigma$, by $p_{div,\Sigma_C}^*$ the critical divergence load of $\Sigma_C$ and by $p_{sw}^*$ the critical value of $p$ for the second order work criterion, we then get

$$p_{sw}^* \le p_{div,C}^* \quad \forall C \quad (4)$$

even though $C = \emptyset$ and the value $p_{sw}^*$ is optimal. Finally, $p_{sw}^* \le p_{div}^*$ (meaning $C = \emptyset$) implying that, as long as $p < p_{sw}^*$ neither the free system $\Sigma$ nor any constrained system $\Sigma_C$ may be divergence unstable.
We are here now concerned by the same issue for flutter (in)stability. We then have to compare the critical flutter load $p_{fl}^*$ of $\Sigma$ with the critical flutter load $p_{fl,C}^*$ of any constrained system $\Sigma_C$. If flutter ki.s.s. is conditional, we have also to find out an eventual value $p_{k,fl}^*$ such that

$$p_{k,fl}^* \le p_{fl,C}^* \quad \forall C \qquad (5)$$

$(p_{fl}^*$ is the critical flutter load of $\Sigma$ and $p_{fl,C}^*$ is the critical flutter load of $\Sigma_C$). Only the case of a free 3 d.o.f system $\Sigma = \Sigma_{free}$ is investigated and, in order to preserve the possibility for the constrained system $\Sigma_C$ to be flutter unstable, only one constraint ($\ell = 1$) is possible.

## 1.2. Flutter ki.s.s. with $n = 3$

Let $u = u(p) \in \mathcal{L}(E)$ be the morphism of the euclidean $n$-dimensional vector space $E = \mathbb{R}^n$ with matrix $\tilde{K}(p)$ in the canonical basis of $E$. Generally, the flutter instability of the system $\Sigma$ means that the operator $u$ fails to be $\mathbb{R}$-diagonalizable.

According to [1], the ki.s.s. is relative to the so-called compressions $u_F$ of $u$ to the subspaces $F$ of $E$ defined by $u_F = p_F \circ u_{|F} \in \mathcal{L}(F)$ where $u_{|F}$ is the restriction of $u$ to $F$ and $p_F$ is the orthogonal projection on $F$. According to the investigated criterion, the ki.s.s. is equivalent to study if the considered property of $u$ is preserved for all its compressions $u_F$.

For divergence stability, the ki.s.s. issue and its solution can be so reformulated: the compressions of an invertible linear map $u = u(p) \in \mathcal{L}(E)$ to any (non nil) subspace $F$ remain invertible if and only if the symmetric part $u_s(p)$ of $u(p)$ is definite. Because for $p = 0$, $u(0)$ is supposed symmetric definite positive, by continuity we are led to the second order work criterion. Recall that the compression of an operator $u$ naturally arises for computing its numerical range $W(u)$. The relationship between the numerical range $W(K)$ of the stiffness matrix $K$ and the second order work is direct: $K$ satisfies the second order work criterion means that $W(K) \subset \mathbb{R}_+^*$. For compressions and numerical range of matrices and operators see for example [9].

Supposed now $u(p) \in \mathbb{R}$-diagonalizable with, to simplify, only simple eigenvalues. The flutter ki.s.s. issue then consists to know if there is a vector subspace $F$ of $E$ such that $u_F = p_F \circ u_{|F} \in \mathcal{L}(F)$ is no more $\mathbb{R}$-diagonalizable and eventually to find such a candidate $F_{k,fl}$. This issue in its great generality is very complicated first because there no general convenient criterion of $\mathbb{R}$ diagonalizability and secondly because, even with a practicable algebraic criterion for low dimensions, the issue remains, as we will see below, a real challenge.

Suppose that $\dim(E) = 3$ (3 d.o.f. system $\Sigma = \Sigma_{free}$) and that the system is constrained by only one kinematic constraint. This constraint is described by a vector $e_3$ that can be chosen on the sphere $S(E)$. In fact, $e_3$ or $-e_3$ represents the same constraint showing that in this case the good geometric structure is the projective space. However, in the general case, the constraint is multidimensional and represented by the vector space $F^\perp$. That means that the well adapted structure to investigate the general issue is, as already mentioned in [1] and in the introduction, the one of Grassmann manifold $Gr_{m,n}(\mathbb{R}) = Gr_m(\mathbb{R}^n) = Gr_m(E)$ of all $m$-dimensional subspaces of $E$ which is a $m(n-m)$ dimensional compact manifold. To conclude this paragraph, let us remind that the $\mathbb{R}$-diagonalizability used criterion for $u_F$ when $\dim(F) = m = 2$ is

$$\Delta(\chi_{u_F}) = \text{tr}^2(u_F) - 4 \det(u_F) > 0 \qquad (6)$$

where $\chi_{u_F}$ is the characteristic polynomial of $u_F$ and that leads to the following flutter ki.s.s criterion:

$$\min_{F \in Gr_2(E)} \text{tr}^2(u_F) - 4 \det(u_F) > 0. \qquad (7)$$

Because of the compactness of $Gr_2(E)$ and the continuity of

$$\Phi : Gr_2(E) \to \mathbb{R} \\ F \mapsto \Phi(F) = \text{tr}^2(u_F) - 4 \det(u_F) \qquad (8)$$
the minimum exists and is reached for an element $F_{k,fl}$. The corresponding constraint is then given by any vector $e_3 \in F_{k,fl}^{\perp} \cap S(E)$.

## 2. First results

### 2.1. Geometric considerations and preliminary calculations

The aim of this paragraph is to transform (8) in order to solve (7). Indeed, the natural way to calculate $F_{k,fl}$ consists on differentiating $\phi$ for getting the critical points. The derivative of $\Phi$ is however difficult to be evaluated because the “points” are vector spaces and we then start by transforming $\Phi$. To do it, we better view $Gr_2(E)$ as the quotient space of the Stiefel manifold $St_2(E)$ of all 2-uplets $(e_1, e_2)$ of orthonormal vectors of $E$ by the orthogonal group $O_2(\mathbb{R})$. We then will describe any point $F$ of $Gr_2(E)$ and any function of $F$ like $\Phi$ by its expression as function of any element $(e_1, e_2)$ of $St_2(E)$ such that the vector space spanned by $(e_1, e_2)$ is $F$ without forgetting that it may be independent of this choice because of the quotient by $O_2(\mathbb{R})$. Remark that, by this way, $St_2(E)$ appears as a principal fiber bundle with $O_2(\mathbb{R})$ as group and with $Gr_2(E)$ as base space: $Gr_2(E) = St_2(E)/O_2(\mathbb{R})$.

Here, $\dim St_2(E) = 3$, $\dim O_2(\mathbb{R}) = 1$ and $\pi: St_2(E) \to Gr_2(E)$ is the projection map of the total space $St_2(E)$ of the principal fiber bundle on its base space $Gr_2(E) = St_2(E)/O_2(\mathbb{R})$ that at each family $(e_1, e_2)$ of two orthonormal vectors of $E$ associates the vector space $\pi(e_1, e_2)$ spanned by these vectors. Reciprocally, if $F \in Gr_2(E) = St_2(E)/O_2(\mathbb{R})$, the set $\pi^{-1}(F)$ is the fiber over $F$ built by all the orthonormal bases $(e_1, e_2)$ of $F$. Let $(e_1, e_2)$ be an orthonormal basis of $F = \pi((e_1, e_2))$. We use the letter $\phi$ instead of the same capital letter $\Phi$ to refer to the function of the variables in $St_2(E)$ so that $\Phi \circ \pi = \phi$. $\phi$ is a lift of $\Phi$. Two expressions of $\Phi$ that will be used are given by:

**Lemma 1.**

$$ \begin{aligned} \Phi(F) = \phi((e_1, e_2)) &= 4(u(e_1) | e_2)(u(e_2) | e_1) + ((u(e_1) | e_1) - (u(e_2) | e_2))^2 \quad (9) \\ &= ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 - 4(u_a(e_1) | e_2)^2. \quad (10) \end{aligned} $$

**Proof.** The following transformations hold:

$$ \mathrm{tr}(u_F) = \sum_{i=1}^{2} (u_F(e_i) | e_i) = \sum_{i=1}^{2} (p \circ u(e_i) | e_i) = \sum_{i=1}^{2} (u(e_i) | e_i) \quad (11) $$

because $p$ is self-adjoint and

$$ \begin{aligned} \det(u_F) &= (u_F(e_1) | e_1)(u_F(e_2) | e_2) - (u_F(e_1) | e_2)(u_F(e_2) | e_1) \\ &= (p \circ u(e_1) | e_1)(p \circ u(e_2) | e_2) - (p \circ u(e_1) | e_2)(p \circ u(e_2) | e_1) \\ &= (u(e_1) | e_1)(u(e_2) | e_2) - (u(e_1) | e_2)(u(e_2) | e_1) \end{aligned} \quad (12) $$

for the same reasons. Straightforward calculations show then that:

$$ \Phi(F) = \phi((e_1, e_2)) = 4(u(e_1) | e_2)(u(e_2) | e_1) + ((u(e_1) | e_1) - (u(e_2) | e_2))^2 $$

which is exactly (9). Here $(e_1, e_2)$ is viewed as an element of $\pi^{-1}(F) \subset St_2(E)$. Direct calculations may show that this expression of $\phi((e_1, e_2))$ does not depend on the choice of $(e_1, e_2)$ as orthonormal basis of $F$ meaning in the fiber over $F$ and then justifying the notation $\Phi(F)$. That may be directly checked at each step of the calculations but we will not do it.

To better understand the flutter as a competition between the symmetric part $u_s$ of $u$ and its skew symmetric part $u_a$ and because the second order work criterion involves $u_s$, we now transform the last expression
by using the symmetry of $u_s$ and the skew symmetry of $u_a$. Calculations give:

$$ \begin{aligned} (u(e_1) | e_2)(u(e_2) | e_1) &= ((u_s(e_1) | e_2) + (u_a(e_1) | e_2))((u_s(e_2) | e_1) + (u_a(e_2) | e_1)) \\ &= (u_s(e_1) | e_2)^2 - (u_a(e_1) | e_2)^2 \end{aligned} $$

leading to:

$$ \begin{aligned} \Phi(F) = \phi((e_1, e_2)) &= 4((u_s(e_1) | e_2)^2 - (u_a(e_1) | e_2)^2) + ((u_s(e_1) | e_1) - (u_s(e_2) | e_2))^2 \\ &= ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 - 4(u_a(e_1) | e_2)^2 \end{aligned} $$

which is exactly (10). ■

It may be checked again that this expression of $\phi((e_1, e_2))$ does not depend on the choice of $(e_1, e_2)$ but only on its equivalence class under the group action by $O_2(\mathbb{R})$ which justify the expression $\Phi(F)$.

But there is a very significant choice. Let be $e_3$ any of the both unit vectors of $F^{\perp}$. Suppose that $e_3 \notin$ ker $u_a$ (this last case will be handled separately and when $n = 3$, it is a one dimensional vector space). Then, the orthogonal space $(F^{\perp})^{\perp} = F$ is spanned by $(u_a(e_3), u_a^2(e_3))$ so that one can choose $e_1 = \frac{u_a(e_3)}{\|u_a(e_3)\|}$, $e_2 = e_3 \wedge \frac{u_a(e_3)}{\|u_a(e_3)\|}$.

Define $\Delta_1(e_3) = \det(e_3, u_a(e_3), u_s(e_3))$, $\Delta_k(e_3) = \det(e_3, u_a(e_3), u_a^k(e_3))$ for $k \ge 2$. Straightforward calculations give:

$$ \begin{aligned} \Phi(F) = \phi((e_1, e_2)) &= h(e_3) = ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 - 4(u_a(e_1) | e_2)^2 \\ &= (\mathrm{Tr} \, u_s - (u_s(e_3) | e_3))^2 - 4 \left( \frac{u_a^2(e_3) | e_3 \wedge u_a(e_3)}{\|u_a(e_3)\|^2} \right)^2 \\ &= (\mathrm{Tr} \, u_s - (u_s(e_3) | e_3))^2 - 4 \left( \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} \right)^2 \end{aligned} \quad (13) $$

4 intrinsic quantities are involved in the issue: the three real eigenvalues of $u_s$: $\alpha_1 \le \alpha_2 \le \alpha_3$ with $(v_i)_i$ an adapted orthonormal basis of $E$ ($u_s(v_i) = \alpha_i v_i$ for all $i$) and $-\beta^2 (\beta > 0)$ the unique not nil eigenvalue of the symmetric linear map $u_a^2$. Indeed, because $n = 3$, rank $u_a = 2$, ker $u_a = \ker u_a^2$ and $(\ker u_a)^\perp = \mathrm{Im}u_a = \mathrm{Im}u_a^2 = E_{-\beta^2} = \langle w_1, w_2 \rangle$ with $w_1$ and $w_2$ two orthonormal eigenvectors of $u_a^2$ so that $u_a^2(w_i) = -\beta^2 w_i$ for $i = 1, 2$. We then deduce, after having chosen $w_3 \in \ker u_a$, that $u_a(w_1) = \beta w_2$, $u_a(w_2) = -\beta w_1$, $u_a(w_3) = 0$ ($w_i)_i$ is still an orthonormal basis of $E$).

The quantities $\Delta_k(e_3)$ will play a significant rule and they can be explicitly evaluated:

**Proposition 1.** 1. For all $x \in E$ and for all $k \ge 1$, $\Delta_{2k+1}(x) = 0$.

2. For all $x \in S(E)$, $\Delta_{2k}(x) = (-1)^{k-1}\beta^2(k-1)\Delta_2(x)$.

3. For all $x \in S(E)$, $\Delta_{2k}(x) = (-1)^{k-1}\beta^2(k-1)\|u_a(x)\|^2 \sqrt{\beta^2 - \|u_a(x)\|^2}$.

**Proof.** These results are obvious if $x \in \ker u_a$ and because $\Delta_k$ is a 3 homogeneous function of $x$, we may suppose that $x \in S(E) \setminus \ker u_a$. Moreover, the minimal polynomial of $u_a^2$ is $\pi_{u_a^2} = X(X + \beta^2)$ meaning that $u_a^4 = -\beta^2 u_a^2$ that leads to the second assertion 2.

Moreover,

$$ \begin{aligned} \Delta_3(x) &= \det(x, u_a(x), u_a^3(x)) = (x \wedge u_a(x) | u_a^3(x)) \\ &= -\|u_a(x)\|^2(x \wedge u_a(x) | u_a(x)) + \frac{\Delta_2(x)}{\|u_a(x)\|^2}(x \wedge u_a(x) | u_a(x \wedge u_a(x))) \\ &= 0 + 0 = 0 \end{aligned} $$

that proves the first assertion.
Finally, the last assertion comes from the following calculation of $\Delta_2(x)$. By using the orthonormal basis $(x, \frac{u_a(x)}{\|u_a(x)\|}, x \wedge \frac{u_a(x)}{\|u_a(x)\|})$, $u_a^2(x)$ reads:

$$ 
\begin{aligned}
u_a^2(x) &= (u_a^2(x) | x)x + \frac{(u_a^2(x) | u_a(x))}{\|u_a(x)\|^2} u_a(x) + \frac{(u_a^2(x) | x \wedge u_a(x))}{\|u_a(x)\|^2} x \wedge u_a(x) \\
&= (u_a^2(x) | x)x + \frac{\Delta_2(x)}{\|u_a(x)\|^2} x \wedge u_a(x) \\
&= -\|u_a(x)\|^2 x + \frac{\Delta_2(x)}{\|u_a(x)\|^2} x \wedge u_a(x).
\end{aligned} 
\quad (14) $$

Calculating now the square of the norm of the right-hand side and the left-hand side of the last equation, we get:

$$ 
\begin{aligned}
\|u_a^2(e_3)\|^2 &= \| - \|u_a(e_3)\|^2 e_3 + \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} e_3 \wedge u_a(e_3) \|^2 \\
(u_a^2(e_3) | u_a^2(e_3)) &= \|u_a(e_3)\|^4 + \left( \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} \right)^2 \|e_3 \wedge u_a(e_3)\|^2 + \\
&\quad - 2\|u_a(e_3)\|^2 \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} (e_3 | e_3 \wedge u_a(e_3)) \\
(u_a^4(e_3) | e_3) &= \|u_a(e_3)\|^4 + \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2} \\
-\beta^2(u_a^2(e_3) | e_3) &= \|u_a(e_3)\|^4 + \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2} \\
\beta^2(u_a(e_3) | u_a(e_3)) &= \|u_a(e_3)\|^4 + \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2} \\
\beta^2\|u_a(e_3)\|^2 &= \|u_a(e_3)\|^4 + \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2}
\end{aligned} 
$$

leading to

$$ \Delta_2^2(e_3) = \|u_a(e_3)\|^4 (\beta^2 - \|u_a(e_3)\|^2) \quad (15) $$

and because $\Delta_2(e_3) \ge 0$:

$$ \Delta_2(e_3) = \|u_a(e_3)\|^2 \sqrt{\beta^2 - \|u_a(e_3)\|^2} \quad (16) $$

which proves the last assertion. ■

## 2.2. Sufficient conditions

First and foremost, remark that in the degenerated case where the second order work criterion (SOWC) fails with a 2 dimensional isotropic cone $C$, choosing $F \subset C$ and the constraint in $F^\perp$ destabilizes the system and flutter ki.s.s. fails. Indeed, we then get $\text{tr}(u_{s,F}) = 0$ and, according to (10), we get $\Phi(F) = -4(u_a(e_1) | e_2)^2 \le 0$.

On the contrary, suppose now to simplify the reasoning, that the SOWC holds meaning here that $\alpha_1 > 0$. Then, without calculating the minimum of $\Phi$ (or $\phi, h$), a sufficient flutter ki.s.s. condition may be got.

**Proposition 2.** As long as

$$ \alpha_1 + \alpha_2 > 2\beta \quad (17) $$
the flutter k.i.s.s. is ensured. Moreover, if for $p = 0$ the system is elastic conservative stable then, the flutter k.i.s.s. is ensured on $[0, p_{skf}]$ where $[p_{skf}]$ is minimal positive root of

$$ \alpha_1(p) + \alpha_2(p) - 2\beta(p) = 0. \qquad (18) $$

**Proof.** Because of well-known results about Rayleigh's quotient for $u_s$, $\alpha_1 \le (u_s(x) | x) \le \alpha_3$ for all unit vector $x$ and the extrema are respectively reached for the eigenvectors $v_1$ (minimum) and $v_3$ (maximum) associated to $\alpha_1$ and $\alpha_3$. Then

$$ \alpha_1 + \alpha_2 \le (u_s(e_1) | e_1) + (u_s(e_2) | e_2) = \text{tr}(u_s) - (u_s(e_3) | e_3) = \sum_{i=1}^{3} \alpha_i - (u_s(e_3) | e_3) \le \alpha_2 + \alpha_3 $$

with a minimum when $e_3 = v_3$ and a maximum when $e_3 = v_1$.

Moreover,

$$ 0 \le |(u_a(e_1) | e_2)| \le \beta $$

with a minimum when $e_3 \in \ker u_a^\perp$ and a maximum when $e_3 = w_3$. Then:

$$ \Phi(F) \ge (\alpha_1 + \alpha_2)^2 - 4\beta^2 \qquad (19) $$

and we deduce that a flutter k.i.s.s. sufficient condition reads as (17) meaning that, as long as the arithmetic mean of the both lowest eigenvalues of $u_s$ is greater than the square root of the not nil eigenvalue of $u_a^2$, no additional kinematic constraint may destabilize the system $\Sigma$ and the flutter k.i.s.s. is ensured. Moreover, suppose as usually that, for $p = 0$, the system is elastic conservative stable. Then, $\alpha_1(0) > 0, \alpha_2(0) > 0$ and $\beta(0) = 0$. Thus, by continuity the minimal positive value $p_{skf}$ root of $\alpha_1(p) + \alpha_2(p) - 2\beta(p) = 0$ is $>0$. On $[0, p_{skf}]$, the flutter k.i.s.s. is ensured. ■

Suppose now that $\alpha_1 + \alpha_2 \le 2\beta$. Because both terms in competition are reached for $e_3 = v_3$ and for $e_3 = w_3$ and because $v_3 \ne w_3$, there is no chance in order that this equality should be realized by a convenient constraint and the flutter k.i.s.s. can be still ensured. The sufficient condition (17) is then not necessary nor optimal and we now tackle the issue of necessary and sufficient flutter k.i.s.s. conditions.

### 3. Calculation of the extrema of $\Phi$ and $h$

The aim of this section is then to calculate the minimum of $\Phi$, $\phi$, $h$. It is a significant challenge because of the nature of the variables, the non convexity of the function and the deep non linearity of the issue. It will be done separately for $\phi$ and $h$, the one for $\Phi$ resulting from those last both. It allows first to validate the results and secondly to highlight the power of the geometric tools like Grassmann or Stiefel manifolds which are the good tools for generalizing up to any dimension $n \ge 4$ the problem supposed here tridimensional. The two conditions defining the critical points are themselves nonlinear because the functions $\Phi$, $\phi$, $h$ are not quadratic. More specifically, they are roughly speaking 4-homogeneous. So, there is no hope of giving the solution by an algorithm involving only linear algebra instructions like for the divergence k.i.s.s and the second order work criterion. The conditions have however a very nice expression and a significant geometrical meaning. Let us start.

#### 3.1. Extremum of $\Phi$

The first step is to find the tangent spaces $T_{(e_1,e_2)}St_2(E)$ and $T_{\pi(e_1,e_2)}Gr_2(E) = T_FGr_2(E)$ for all $(e_1, e_2) \in St_2(E)$ and $\pi(e_1, e_2) = F \in Gr_2(E)$.
Usual results about Stiefel manifolds show that

**Lemma 2.** For any orthonormal family $(e_1, e_2) \in St_2(E)$,

$$T_{(e_1, e_2)} St_2(E) = \{(\epsilon_1 = se_2 + t_1 e_3, \epsilon_2 = -se_2 + t_2 e_3) \mid s, t_1, t_2 \in \mathbb{R}\}$$

and that

$$T_{\pi(e_1, e_2)} Gr_2(E) = T_F Gr_2(E) = \pi^T(T_{(e_1, e_2)} St_2(E)) = \{(\epsilon_1 = t_1 e_3, \epsilon_2 = t_2 e_3) \mid t_1, t_2 \in \mathbb{R}\}$$

where $e_3$ is any unit vector orthogonal to $F = \pi((e_1, e_2))$.

With the same notations, critical points of $\Phi$ are characterized by the following.

**Proposition 3.** $F = \pi((e_1, e_2))$ is a critical point of $\Phi$ if both following relations hold if

$$(u_s(e_3) \mid u_a(e_3)) = 0 \qquad (20)$$

$$(\operatorname{tr} u_s - (e_3 \mid u_s(e_3)))\Delta_1(e_3) + 2\Delta_2(e_3) = 0 \qquad (21)$$

for any normalized vector $e_3$ such that $F = \langle e_3 \rangle^{\perp}$.

**Proof.** Usual derivative calculations give:

$$\phi'(e_1, e_2)(\epsilon_1, \epsilon_2) = 2((u_s(e_1) \mid e_1) + (u_s(e_2) \mid e_2))(2(u_s(e_1) \mid \epsilon_1) + 2(u_s(e_2) \mid \epsilon_2)) - 8(u_a(e_1) \mid e_2)((u_a(\epsilon_1) \mid e_2) + (u_a(e_1) \mid \epsilon_2)) \quad \text{and thus}$$

$$\Phi'(\pi((e_1, e_2))(t_1 e_3, t_2 e_3)) = 4((u_s(e_1) \mid e_1) + (u_s(e_2) \mid e_2))(t_1(u_s(e_1) \mid e_3) + t_2(u_s(e_2) \mid e_3)) - 8(u_a(e_1) \mid e_2)(t_1(u_a(e_3) \mid e_2) + t_2(u_a(e_1) \mid e_3))$$

$F = \pi((e_1, e_2))$ is a critical point of $\Phi$ if and only if $\Phi'(\pi((e_1, e_2))(t_1 e_3, t_2 e_3)) = 0$ for all $t_1, t_2 \in \mathbb{R}$. For successively $(t_1, t_2) = (1, 0)$ and $(t_1, t_2) = (0, 1)$, it leads to:

$$((u_s(e_1) \mid e_1) + (u_s(e_2) \mid e_2))(u_s(e_1) \mid e_3) = 2(u_a(e_1) \mid e_2)(u_a(e_3) \mid e_2) \qquad (22)$$

$$((u_s(e_1) \mid e_1) + (u_s(e_2) \mid e_2))(u_s(e_2) \mid e_3) = 2(u_a(e_1) \mid e_2)(u_a(e_1) \mid e_3) \qquad (23)$$

(the reader may check that these relations are really invariant by any rotation about $e_3 \in F^{\perp}$). These relations are the both (nonlinear) conditions on $(e_1, e_2)$ but actually on $F$ defining the critical points of $\Phi$. They may however be transformed by using a unit vector $e_3 \in F^{\perp}$.

From $u_s(e_i) = \sum_{k=1}^3 (u_s(e_i) \mid e_k)e_k$ and $u_a(e_i) = \sum_{k=1, k \neq i}^3 (u_a(e_i) \mid e_k)e_k$, we deduce

$$(u_s(e_1) \mid u_a(e_1)) = (u_s(e_1) \mid e_2)(u_a(e_1) \mid e_2) + (u_s(e_1) \mid e_3)(u_a(e_1) \mid e_3) \qquad (24)$$

$$(u_s(e_2) \mid u_a(e_2)) = (u_s(e_2) \mid e_1)(u_a(e_2) \mid e_1) + (u_s(e_2) \mid e_3)(u_a(e_2) \mid e_3). \qquad (25)$$

But from (22) and (23), we get

$$\frac{(u_s(e_1) \mid e_3)}{(u_s(e_2) \mid e_3)} = \frac{(u_a(e_3) \mid e_2)}{(u_a(e_1) \mid e_3)} \left( = \frac{2(u_a(e_1) \mid e_2)}{(u_s(e_1) \mid e_1) + (u_s(e_2) \mid e_2)} \right)$$

or

$$(u_s(e_1) \mid e_3)(u_a(e_1) \mid e_3) = (u_a(e_3) \mid e_2)(u_s(e_2) \mid e_3)$$

or still

$$(u_s(e_1) \mid e_3)(u_a(e_1) \mid e_3) + (u_a(e_2) \mid e_3)(u_s(e_2) \mid e_3) = 0.$$

But as $u_a$ is skew symmetric and $u_s$ symmetric, from (24) and (25)

$$(u_s(e_1) \mid u_a(e_1)) + (u_s(e_2) \mid u_a(e_2)) = (u_s(e_1) \mid e_3)(u_a(e_1) \mid e_3) + (u_s(e_2) \mid e_3)(u_a(e_2) \mid e_3).$$
Thus

$$ (u_s(e_1) | u_a(e_1)) + (u_s(e_2) | u_a(e_2)) = 0. $$

But

$$ \begin{aligned} (u_s(e_3) | u_a(e_3)) &= (u_s(e_3) | e_1)(u_a(e_3) | e_1) + (u_s(e_3) | e_2)(u_a(e_3) | e_2) \\ &= -(u_s(e_1) | e_3)(u_a(e_1) | e_3) - (u_s(e_2) | e_3)(u_a(e_2) | e_3). \end{aligned} $$

Thus

$$ (u_s(e_3) | u_a(e_3)) = 0 \qquad (26) $$

which is exactly (20).

Because one may choose any unit vector orthogonal $e_3$ for $e_1$, one takes $e_1 = \frac{u_a(e_3)}{\|u_a(e_3)\|}$ and $e_2 = e_3 \wedge e_1 = \frac{e_3 \wedge u_a(e_3)}{\|u_a(e_3)\|}$. (22) is then reduced to $0 = 0$ and calculations give successively:

$$ \begin{aligned} (u_s(e_2) | e_3) &= (e_2 | u_s(e_3)) = \frac{\det(e_3, u_a(e_3), u_s(e_3))}{\|u_a(e_3)\|} = \frac{\Delta_1(e_3)}{\|u_a(e_3)\|} \\ (u_a(e_1) | e_2) &= \left( \frac{u_a^2(e_3)}{\|u_a(e_3)\|} \middle| \frac{e_3 \wedge u_a(e_3)}{\|u_a(e_3)\|} \right) = \frac{\det(e_3, u_a(e_3), u_a^2(e_3))}{\|u_a(e_3)\|^2} = \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} \end{aligned} \qquad (27) $$

and

$$ (u_a(e_1) | e_3) = -(e_1 | u_a(e_3)) = - \left( \frac{u_a(e_3)}{\|u_a(e_3)\|} \middle| u_a(e_3) \right) = - \|u_a(e_3)\|. $$

Thus

$$ (u_a(e_1) | e_2)(u_a(e_1) | e_3) = - \frac{\det(e_3, u_a(e_3), u_a^2(e_3)) \|u_a(e_3)\|}{\|u_a(e_3)\|^2} = - \frac{\det(e_3, u_a(e_3), u_a^2(e_3))}{\|u_a(e_3)\|} = - \frac{\Delta_2(e_3)}{\|u_a(e_3)\|} $$

and finally

$$ \begin{aligned} (23) \Leftrightarrow & ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))\Delta_1(e_3) + 2\Delta_2(e_3) = 0 \quad \text{or} \\ & (\text{tr } u_s - (e_3 | u_s(e_3)))\Delta_1(e_3) + 2\Delta_2(e_3) = 0 \end{aligned} \qquad (28) $$

the last equation being exactly (21). ■

By evaluating $\Phi$ at a critical point, we get the following flutter ki.s.s condition:

**Proposition 4.** Flutter ki.s.s holds as long as $\frac{\Delta_1^2(e_3)}{\|u_a(e_3)\|^4} < 1$ or when

$$ | \Delta_1(e_3) = \det(e_3, u_a(e_3), u_s(e_3)) | < \|u_a(e_3)\|^2 = \det(e_3, u_a(e_3), e_3 \wedge u_a(e_3)) \qquad (29) $$

or when

$$ \left| \det \left( e_3, \frac{u_a(e_3)}{\|u_a(e_3)\|}, \frac{u_s(e_3)}{\|u_s(e_3)\|} \right) \right| < \frac{\|u_a(e_3)\|}{\|u_s(e_3)\|} \qquad (30) $$

for all units $e_3$ such that $F = \pi((e_1, e_2)) = \langle e_3 \rangle^\perp$ is a critical point of $\Phi$.

Geometrically speaking, (32) means that the volume of the parallelepiped built on $(e_3, u_a(e_3), u_s(e_3))$ is lower than that built on $(e_3, u_a(e_3), e_3 \wedge u_a(e_3))$ (equal to $\|u_a(e_3)\|^2$).

**Proof.** Remark that from (27) and (28), we deduce that, at a critical point $F = \pi((e_1, e_2)) = \langle e_3 \rangle^\perp$, $\Phi$ reads:
$$
\begin{align}
\Phi \circ \pi((e_1, e_2)) &= ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 - 4(u_a(e_1) | e_2)^2 \nonumber \\
&= ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 - \left( \frac{2\Delta_2(e_3)}{\|u_a(e_3)\|^2} \right)^2 \nonumber \\
&= ((u_s(e_1) | e_1) + (u_s(e_2) | e_2))^2 \left( 1 - \frac{\Delta_1^2(e_3)}{\|u_a(e_3)\|^4} \right). \tag{31}
\end{align}
$$

Flutter ki.s.s holds as long as $\frac{\Delta_1^2(e_3)}{\|u_a(e_3)\|^4} < 1$ or when

$$
| \Delta_1(e_3) = \det(e_3, u_a(e_3), u_s(e_3)) | < \|u_a(e_3)\|^2 = \det(e_3, u_a(e_3), e_3 \wedge u_a(e_3)) \quad (32)
$$

that may be rewritten as

$$
\left| \det \left( e_3, \frac{u_a(e_3)}{\|u_a(e_3)\|}, \frac{u_s(e_3)}{\|u_s(e_3)\|} \right) \right| < \frac{\|u_a(e_3)\|}{\|u_s(e_3)\|} \quad (33)
$$

for all unit $e_3$ such that $F = \pi((e_1, e_2)) = \langle e_3 \rangle^{\perp}$ is a critical point of $\Phi$ namely solutions of (20) and (28)
(or (21)) or equivalently of (22) and (23). ■

3.2. Extremum of h

To validate the previous results and to present the calculations with more usual tools of differential ge-
ometry, we use the parametrization of the problem by the sphere S(E) of unit vectors of E meaning by the
function h defined by (13) that is now recalled:

$$
h(e_3) = (\operatorname{tr} u_s - (u_s(e_3) | e_3))^2 - 4 \left( \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} \right)^2 . \quad (34)
$$

The aim of this paragraph is then to calculate the critical points of $h$ on the sphere $S(E)$.

Put $g(e_3) = \operatorname{tr} u_s - (u_s(e_3) | e_3)$ and $f(e_3) = \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2}$ so that $h(e_3) = g^2(e_3) - 4f^2(e_3)$. Straightforward calculations give:

$$
\begin{align}
h'(e_3)(\epsilon) &= 2g(e_3)g'(e_3)(\epsilon) - 8f(e_3)f'(e_3)(\epsilon) \tag{35} \\
g'(e_3)(\epsilon) &= -2(e_3 | u_s(\epsilon)) \nonumber \\
f'(e_3)(\epsilon) &= \frac{\Delta_2'(e_3)(\epsilon)\|u_a(e_3)\|^2 - 2\Delta_2(e_3)(u_a(e_3) | u_a(\epsilon))}{\|u_a(e_3)\|^4} \nonumber \\
\Delta_2'(e_3)(\epsilon) &= \det(\epsilon, u_a(e_3), u_a^2(e_3)) + \det(e_3, u_a(\epsilon), u_a^2(e_3)) + \det(e_3, u_a(e_3), u_a^2(\epsilon)) \nonumber
\end{align}
$$

for all $\epsilon \in T_{e_3}S(E) = \langle e_3 \rangle^{\perp}$ and for getting the critical points $e_3$, it is necessary and sufficient to verify
$h'(e_3)(\epsilon) = 0$ on a basis of $T_{e_3}S_3(E) = e_3^{\perp}$ meaning for $\epsilon = u_a(e_3)$ and $\epsilon = e_3 \wedge u_a(e_3)$ ($e_3 \not\in \ker u_a$). We do
it in the following two subsections.

3.2.1. $\epsilon = u_a(e_3)$

**Proposition 5.** $h'(e_3)(u_a(e_3)) = 0$ *implies relation (20).*

**Proof.** One successively finds:

$$
\begin{align*}
g'(e_3)(\epsilon) &= g'(e_3)(u_a(e_3)) \\
&= -2(e_3 | u_s(u_a(e_3))) \\
&= -2(u_s(e_3) | u_a(e_3))
\end{align*}
$$
$$
\begin{align*}
\Delta'_2(e_3)(\epsilon) &= \Delta'_2(e_3)(u_a(e_3)) \\
&= \det(u_a(e_3), u_a(e_3), u_a^2(e_3)) + \det(e_3, u_a^2(e_3), u_a^2(e_3)) + \det(e_3, u_a(e_3), u_a^3(e_3)) \\
&= \det(e_3, u_a(e_3), u_a^3(e_3)) \\
&= \Delta_3(e_3) = 0
\end{align*}
$$

because of Proposition 1. We deduce that:

$$
\begin{align*}
f'(e_3)(\epsilon) &= f'(e_3)(u_a(e_3)) \\
&= \frac{\Delta'_2(e_3)(u_a(e_3))\|u_a(e_3)\|^2 - 2\Delta_2(e_3)(u_a(e_3) \mid u_a^2(e_3))}{\|u_a(e_3)\|^4} \\
&= \frac{-2\Delta_2(e_3)(u_a(e_3) \mid u_a^2(e_3))}{\|u_a(e_3)\|^4} = 0
\end{align*}
$$

because of the skew-symmetry of $u_a$. From (35), we then deduce $g'(e_3) = 0$ ($g(e_3) \neq 0$) or $(u_s(e_3) | u_a(e_3)) = 0$. We find again (20). ■

$$
3.2.2. \quad \epsilon = e_3 \wedge u_a(e_3)
$$

**Proposition 6.** *h'(*e*<sub>3</sub>(*e*<sub>3</sub> ∧ *u*<sub>*a*</sub>(*e*<sub>3</sub>) = 0 implies relation (21).*

**Proof.** The calculations are a little more complicated. Calculations give:

$$
\begin{align*}
g'(e_3)(\epsilon) &= g'(e_3)(e_3 \wedge u_a(e_3)) \\
&= -2(e_3 \mid u_s(e_3 \wedge u_a(e_3))) \\
&= -2(u_s(e_3) \mid e_3 \wedge u_a(e_3)) \\
&= -2\Delta_1(e_3) \\
\Delta'_2(e_3)(\epsilon) &= \Delta'_2(e_3)(e_3 \wedge u_a(e_3)) \\
&= \det(e_3 \wedge u_a(e_3), u_a(e_3), u_a^2(e_3)) + \det(e_3, u_a(e_3 \wedge u_a(e_3)), u_a^2(e_3)) \\
&\quad + \det(e_3, u_a(e_3), u_a^2(e_3 \wedge u_a(e_3))).
\end{align*}
$$

But from (14)

$$
\begin{align*}
\det(e_3 \wedge u_a(e_3), u_a(e_3), u_a^2(e_3)) &= \det \left( e_3 \wedge u_a(e_3), u_a(e_3), -\|u_a(e_3)\|^2 e_3 + \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} e_3 \wedge u_a(e_3) \right) \\
&= -\|u_a(e_3)\|^2 \det(e_3 \wedge u_a(e_3), u_a(e_3), e_3) \\
&= \|u_a(e_3)\|^2 \det(e_3, u_a(e_3), e_3 \wedge u_a(e_3)) \\
&= \|u_a(e_3)\|^4 \tag{36}
\end{align*}
$$

and one more time from (14) and the double cross product formula

$$
\begin{align*}
\det(e_3, u_a(e_3 \wedge u_a(e_3)), u_a^2(e_3)) &= -(e_3 \wedge u_a^2(e_3) \mid u_a(e_3 \wedge u_a(e_3))) \\
&= - \left( e_3 \wedge \left( \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} e_3 \wedge u_a(e_3) \right) \mid u_a(e_3 \wedge u_a(e_3)) \right) \\
&= \left( \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} u_a(e_3) \mid u_a(e_3 \wedge u_a(e_3)) \right) \\
&= - \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} (u_a^2(e_3) \mid e_3 \wedge u_a(e_3)) \\
&= - \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2}. \tag{37}
\end{align*}
$$
From (14), one gets

$$ \frac{\Delta_2(e_3)}{\|u_a(e_3)\|^2} e_3 \wedge u_a(e_3) = \|u_a(e_3)\|^2 e_3 + u_a^2(e_3) $$

and thus successively:

$$ e_3 \wedge u_a(e_3) = \frac{\|u_a(e_3)\|^4}{\Delta_2(e_3)} e_3 + \frac{\|u_a(e_3)\|^2}{\Delta_2(e_3)} u_a^2(e_3) \quad \text{puis} $$

$$ u_a(e_3 \wedge u_a(e_3)) = \frac{\|u_a(e_3)\|^4}{\Delta_2(e_3)} u_a(e_3) + \frac{\|u_a(e_3)\|^2}{\Delta_2(e_3)} u_a^3(e_3) $$

$$ u_a^2(e_3 \wedge u_a(e_3)) = \frac{\|u_a(e_3)\|^4}{\Delta_2(e_3)} u_a^2(e_3) + \frac{\|u_a(e_3)\|^2}{\Delta_2(e_3)} u_a^4(e_3) $$

and

$$ \det(e_3, u_a(e_3), u_a^2(e_3 \wedge u_a(e_3))) = \|u_a(e_3)\|^4 + \frac{\|u_a(e_3)\|^2 \Delta_4(e_3)}{\Delta_2(e_3)} \quad (38) $$

and from Proposition 1

$$ \det(e_3, u_a(e_3), u_a^2(e_3 \wedge u_a(e_3))) = \|u_a(e_3)\|^4 - \beta^2 \|u_a(e_3)\|^2. \quad (39) $$

Finally

$$ \Delta_2'(e_3)(e_3 \wedge u_a(e_3)) = 2\|u_a(e_3)\|^4 - \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^2} - \beta^2 \|u_a(e_3)\|^2. $$

We deduce

$$ \begin{align*} f'(e_3)(e_3 \wedge u_a(e_3)) &= \frac{\Delta_2'(e_3)(e_3 \wedge u_a(e_3)) \|u_a(e_3)\|^2 - 2\Delta_2(e_3)(u_a(e_3) \mid u_a(e_3 \wedge u_a(e_3)))}{\|u_a(e_3)\|^4} \\ &= \frac{\Delta_2'(e_3)(e_3 \wedge u_a(e_3)) \|u_a(e_3)\|^2 + 2\Delta_2(e_3)(u_a^2(e_3) \mid e_3 \wedge u_a(e_3))}{\|u_a(e_3)\|^4} \\ &= \frac{\Delta_2'(e_3)(e_3 \wedge u_a(e_3)) \|u_a(e_3)\|^2 + 2\Delta_2^2(e_3)}{\|u_a(e_3)\|^4} \\ &= 2\|u_a(e_3)\|^2 - \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^4} - \beta^2 + \frac{2\Delta_2^2(e_3)}{\|u_a(e_3)\|^4} \\ &= \frac{\Delta_2^2(e_3)}{\|u_a(e_3)\|^4} + 2\|u_a(e_3)\|^2 - \beta^2 \\ &= \beta^2 - \|u_a(e_3)\|^2 + 2\|u_a(e_3)\|^2 - \beta^2 = \|u_a(e_3)\|^2 \quad \text{again from Proposition 1} \end{align*} $$

$h'(e_3)(e_3 \wedge u_a(e_3)) = 0$ then leads to:

$$ \begin{align*} 0 = h'(e_3)(e_3 \wedge u_a(e_3)) &= 2g(e_3)g'(e_3)(e_3 \wedge u_a(e_3)) - 8f(e_3)f'(e_3)(e_3 \wedge u_a(e_3)) \\ &= -4g(e_3)\Delta_1(e_3) - 8f(e_3)\|u_a(e_3)\|^2 \end{align*} $$
and

$$g(e_3)\Delta_1(e_3) + 2f(e_3)\|u_a(e_3)\|^2 = g(e_3)\Delta_1(e_3) + 2\Delta_2(e_3) = 0 \quad (40)$$

which is exactly (21). ■

### 3.3. Summary of the results

Because the examples will be handled through the parametrization by unit vectors $e_3$, we briefly summarize the results. The critical points $e_{3,c}$ are solutions of the following system:

$$\left\{
\begin{array}{l}
\|e_3\|^2 - 1 = 0 \\
(u_s(e_3) | u_a(e_3)) = 0 \\
(\mathrm{Tr} \, u_s - (u_s(e_3) | e_3))\Delta_1(e_3) + 2\Delta_2(e_3) = 0.
\end{array}
\right.
\quad (41)$$

It is the intersection of three hypersurfaces and it is then built by $p \ge 2$ isolated points $e_{3,c,k}$ in $E$ for $k = 1, \dots, p$. (at least one minimum and one maximum thanks to the compactness of the domain and the continuity of the map). Among these points, there is the wanted absolute minimum $e_{3,m}$ and the flutter ki.s.s. condition then reads:

$$|\Delta_1(e_{3,m}) = \det(e_{3,m}, u_a(e_{3,m}), u_s(e_{3,m}))| < \|u_a(e_{3,m})\|^2 = \det(e_{3,m}, u_a(e_{3,m}), e_{3,m} \wedge u_a(e_{3,m})) \quad (42)$$

and note that $\Delta_1(x) = 2\det(x, u(x), u^*(x))$ for all $x \in E$ which allows us to express the flutter ki.s.s. condition as:

$$|\Delta_1(e_{3,m}) = \det(e_{3,m}, u(e_{3,m}), u^*(e_{3,m}))| < 2\|u_a(e_{3,m})\|^2 = 2\det(e_{3,m}, u_a(e_{3,m}), e_{3,m} \wedge u_a(e_{3,m})). \quad (43)$$

### 3.4. Mechanical consequences

The exact analytic conditions ensuring the flutter ki.s.s. of a flutter stable free system $\Sigma$ are then obtained ((43) with $e_{3,m}$ satisfying (41)) and lead to a critical flutter ki.s.s value $p^*_{k,fl}$ of the load parameter whereas the sufficient condition (18) leads to a value $p^*_{s,k,fl}$ with obviously $p^*_{s,k,fl} \ge p^*_{k,fl}$. We could then say that the flutter ki.s.s is also conditional as for non conservative divergence ki.s.s. But, contrary to the case of divergence ki.s.s. for non conservative systems, the above flutter ki.s.s. conditions are independent of the flutter critical value $p^*_{fl}$ of the free system $\Sigma$ and the both cases $p_{k,fl} < p^*_{fl}$ and $p^*_{k,fl} > p^*_{fl}$ may occur: (5) holds only for all $C$ with $\ell = 1$ but not for $\ell = 0$ (no constraint).

As it has been recalled in Section 1, for a conservative system $\Sigma$ with a divergence critical load $p^*_{div}$, the ki.s.s. is ensured for $p < p^*_{div}$ (universal ki.s.s.) whereas for a non conservative system with a divergence critical load $p^*_{div}$, the ki.s.s. is ensured for $p < p^*_{sw}$ (conditional ki.s.s.) where $p^*_{sw}$ is the second order work critical load that verifies $p^*_{sw} < p^*_{d}$ because $\det K_s(p) \le \det K(p)$.

So, for both conservative and non conservative systems, the divergence ki.s.s. subdomain $[0, p^*_{k,div}]$ (with $p^*_{k,div} = p^*_{div}$ or $p^*_{sw}$ according to the conservative or non conservative cases) is always included in the divergence stability domain $[0, p^*_{div}]$. As a consequence, the divergence stability of every constrained system $\Sigma_C$ also ensures the divergence stability of the free system $\Sigma = \Sigma_{free}$: we may choose $C = \emptyset$ or $\ell = 0$. (In our reasonings, remember that the load is supposed monotonically increasing from $p = 0$ and often implicitly conservative stable for $p = 0$).

On the contrary, the flutter stability of every constrained system $\Sigma_C$ does not more ensure the flutter stability of the free system $\Sigma$: for flutter ki.s.s., to be sure that neither the free system $\Sigma_{free}$ nor any constrained subsystem $\Sigma_C$ may be destabilized, we then have to put $p < \min\{p^*_{fl}, p^*_{k,fl}\}$ without knowing a priori which is the minimum value. Lastly, the value $p^*_{k,fl}$ of flutter ki.s.s. depends on the mass matrix as the value of the critical flutter load parameter $p^*_{fl}$. The following examples will illustrate all these results.
<figure><img src="image_2.png" /><figcaption>Fig. 1. 3 d.o.f. Ziegler system.</figcaption></figure>

# 4. Illustrative examples

## 4.1. The mechanical system: Ziegler's column

We now apply the above results to the usual three d.o.f. Ziegler column $\Sigma$ as the one used in [3] or [1] for investigating the divergence ki.s.s. $\Sigma$ then consists in the three degree of freedom Ziegler system $\Sigma$ as in Fig. 1 made up of three bars OA, AB, BC with $OA = AB = AC = \ell$ linked by three elastic springs of the same stiffness $k$. The nonconservative external action (the circulatory force) is the follower force $\vec{P}$.

The elastic energy of the springs is

$$U = \frac{k}{2}(\theta_1^2 + (\theta_1 - \theta_2)^2 + (\theta_2 - \theta_3)^2) = \frac{k}{2}(2\theta_1^2 + 2\theta_1^2 + \theta_3^2 - 2\theta_1\theta_2 - 2\theta_2\theta_3)$$

and the virtual power of $\vec{P}$ in any configuration $\theta = (\theta_1, \theta_2, \theta_3)$ reads ($P > 0$ in compression):

$$\mathcal{P}_P^* = P\ell(\sin(\theta_3 - \theta_1)\theta_1^* + \sin(\theta_3 - \theta_2)\theta_2^*).$$

Put $p = \frac{P\ell}{k}$ as dimensionless loading parameter and noting that $(0,0,0)$ is the unique equilibrium configuration, the stiffness matrix then reads:

$$K(p) = \begin{bmatrix} 2-p & -1 & p \\ -1 & 2-p & -1+p \\ 0 & -1 & 1 \end{bmatrix}.$$

Two cases of mass matrix will be investigated but we do not systematically give the corresponding mass distribution. Only the second one is associated to a uniform mass distribution. In each case, we give the numerical approximation of the solution $e_3$ (called here $X$ with $X^T = (x_1 x_2 x_3)$) of Eqs. (43) and only for the first case, we give the expanded expressions of the quantities involved in Eqs. (43). We also validate the results by
a direct numerical solution of the initial minimization problem. Because of the non convexity of the criterium function, a specific algorithm is devised. It relies on the conjunction of two classical numerical algorithms. First, a Nelder-Mead downhill simplex technique [10] is used so as to find the minimum of a function $f(X, p)$ with $X$ on the unit sphere for $p$ fixed. Then, this first optimization algorithm is piloted by a dichotomy procedure on $p$ which will converge to the value $p^*$ such that $\min_X f(X, p^*) = 0$. The sphere is parametrized by the usual spherical coordinates $(x_1 = \cos(\psi_1) \cdot \cos(\psi_2), x_2 = \sin(\psi_1) \cdot \cos(\psi_2)$ and $x_3 = \sin(\psi_2)$ with $-\pi \le \psi_1 \le \pi$ and $-\pi/2 \le \psi_2 \le \pi/2$) in order to lead to a minimization problem without constraint.

The critical values $p_{fl}^*, p_{s,k,fl}^*$ and $p_{k,fl}^*$ are then calculated in order to illustrate the above mechanical discussion. For the corresponding value of $X_{k,fl}^T = (a \ b \ c)$ to $p_{k,fl}^*$, the destabilizing kinematic constraint then reads $a\theta_1 + b\theta_2 + c\theta_3 = 0$. The above analytic results give only the first order equations for the critical points and then allows to find a ki.s.s. critical value of the load parameter. In order to have the exact set of equations and inequalities for the minimums points, a second order set of inequalities should be added to select among the critical points the minimum points. But, because of the non convexity of the problem, we have to evaluate $h$ on the critical points or, in the best, on the set of the minimums to select the absolute one. The used simplex technique validates the analytic approach by using a complete different way which avoids any gradient method.

4.2. $M = I_3$

In this first case, $\tilde{K}(p) = K(p)$ and

$$ 
\begin{aligned}
K_a(p) &= \begin{bmatrix} 0 & 0 & 1/2p \\ 0 & 0 & 1/2p \\ -1/2p & -1/2p & 0 \end{bmatrix} \\ 
K_s(p) &= \begin{bmatrix} 2-p & -1 & 1/2p \\ -1 & 2-p & -1+1/2p \\ 1/2p & -1+1/2p & 1 \end{bmatrix}.
\end{aligned}
$$

The vector $e_3$ parametrizing the problem is then the column vector $X$ of the coordinates of $e_3$ in the initial basis.

$$ X = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}. $$

The above quantities then read:

$$ \Delta_1 (X) = \frac{1}{2} x_3 p^2 x_1^2 - x_1 x_2 x_3 p - x_1 p x_3^2 - \frac{1}{2} p x_1^3 - \frac{1}{2} p x_1^2 x_2 - \frac{1}{2} x_1^2 x_3 p + \frac{1}{2} x_1 p x_2^2 + \frac{1}{2} p x_2^3 \\ - \frac{1}{2} x_3 p^2 x_2^2 + \frac{1}{2} x_2^2 x_3 p + x_2 p x_3^2 - \frac{1}{2} x_2 x_3^2 p^2 - \frac{1}{2} p x_3^3 + \frac{1}{2} x_1 x_3^2 p^2 $$

$$ \Delta_2 (X) = -\frac{1}{4} x_1 p^3 x_3^2 - \frac{1}{8} p^3 x_1^3 - \frac{1}{8} p^3 x_1^2 x_2 + \frac{1}{8} x_1 p^3 x_2^2 + \frac{1}{8} p^3 x_2^3 + \frac{1}{4} x_2 p^3 x_3^2 $$

and $h$, that must to be minimized on the sphere, reads:

$$ 
\begin{aligned}
h(X) = & 25 - 20p - 5p^2 x_1^2 - 5p^2 x_2^2 - 4x_1 x_2 x_3^2 + 4x_1 x_2^3 p - 2p^2 x_2^3 x_3 - 6px_2^2 x_3^2 + 8x_2^3 x_3 p + 2x_1 x_3^3 p \\ 
& + x_1^2 x_3^2 p^2 + 2x_2 x_3^3 p + x_2^2 x_3^2 p^2 + 4p^2 + p^2 x_1^4 - 4x_2^4 p + p^2 x_2^4 - 8x_1^3 x_2 + 12x_1^2 x_2^2 \\ 
& + 2p^2 x_1^2 x_2^2 - 2px_1^2 x_3^2 + 4px_1^3 x_2 - 2p^2 x_1^3 x_3 - 8x_1^2 x_2 x_3 - 8x_1^2 p x_2^2 + 4x_1^3 x_3 p \\ 
& + 4x_2 x_3 p^2 + 4x_1 x_3 p^2 - 8px_1 x_2 - 2x_1 x_3 p^2 x_2^2 + 18px_1^2 + 20x_2 x_3 + 20x_1 x_2 - 8x_1 x_2^3 + 4x_2^4 \\ 
& + 4x_1^4 - 2p^2 x_1^2 x_2 x_3 + 4x_1^2 x_2 x_3 p + 2x_1 x_3^2 p^2 x_2 - 4x_1 x_3^2 p x_2 + 4px_3^2 - 4x_1^4 p \\ 
& - 20x_1^2 - 20x_2^2 - 10x_3^2 + 18px_2^2 - 8x_2^3 x_3 + 4x_1^2 x_3^2 - 4x_2 x_3^3 + 8x_2^2 x_3^2 + 8x_1 x_2^2 x_3 \\ 
& - 18x_2 x_3 p - 10x_1 x_3 p + 2p^2 x_1 x_2 + x_3^4.
\end{aligned}
$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) Minimum values of $f(X,p)$ for $X$ on the unit-sphere and $p \in [0, 2]$, (b) 2D map in spherical coordinates $(\psi_1, \psi_2)$ of $f(X,p)$ on the unit-sphere for $p = p_{k,fl}$.</figcaption></figure>

The SOWC is valid as long as $p < 1$ ($p^*sw = 1$) and the flutter appears for the 3 d.o.f (unconstrained) system $\Sigma$ for $p_{fl}^* = 0.515564668$. The necessary and sufficient condition given by (17) leads here to $p_{s,k,fl} = 0.6985965580$. Because $p_{s,k,fl} > p_f$, no kinematic constraint may destabilize by flutter the flutter stable unconstrained system $\Sigma$. The exact value of $p_{k,fl}^*$ is obtained by minimization of $h$ or by solving directly (41) as well and leads to the same value $p_{k,fl}^* = 0.7266116$

$$ (\psi_1 = -0.809337 \text{ and } \psi_2 = 0.332388) $$

fortunately higher than $p_{skf}$. In this case, the flutter ki.s.s. does not imply the flutter stability of the free system $\Sigma = \Sigma_{free}$ and the free system becomes unstable “before” finding a destabilizing kinematic constraint (see Fig. 2).

### 4.3. Uniform mass distribution

No expansion of the quantities like $h$, $\Delta_1$, $\Delta_2$ but only the result values are now given.

For this mass distribution, the mass matrix reads:

$$ M = \begin{pmatrix} \frac{7}{3} & \frac{3}{2} & \frac{1}{2} \\ \frac{3}{2} & \frac{4}{3} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{3} \end{pmatrix}. $$

In this case

$$ \tilde{K}(p) = \begin{pmatrix} 4.786308217 & -1.602849583 p & -7.160390241 + 1.791273952 p & 4.013176468 + 0.233101685 p \\ -7.160390241 + 2.087061729 p & 13.58105314 - 4.615655719 p & -10.76259839 + 4.210033295 p \\ 4.013176468 - 0.6299566260 p & -10.76259839 + 2.040678218 p & 11.86340788 - 2.550725465 p \end{pmatrix}. $$

The SOWC is valid as long as $p < 1$ ($p^*sw = 1$) and the flutter appears for the 3 d.o.f. (unconstrained) system $\Sigma$ for $p_{fl}^* = 1.483549098926$. The necessary and sufficient condition given by (17) leads here to $p_{s,k,fl}^* = 0.9498607071$. Because $p_{s,k,fl}^* < p_{fl}^*$, there is a kinematic constraint that destabilizes by flutter the flutter stable unconstrained system $\Sigma_{free}$. The exact value of $p_{k,fl}^*$ is obtained by minimization of $h$ ($\psi_1 = -0.973998$ and $\psi_2 = 0.666310$) or by solving directly (41) as well and leads to the same value $p_{k,fl}^* = 1.297968$ fortunately again higher than $p_{s,k,fl}^*$. In this case, the flutter ki.s.s. ensures the flutter stability of the free system $\Sigma = \Sigma_{free}$ and the corresponding destabilizing constraint reads (see Fig. 3):

$$ 0.441790153915307\theta_1 - 0.650220378830593\theta_2 + 0.618089737549398\theta_3 = 0. $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Minimum values of $f(X,p)$ for $X$ on the unit-sphere and $p \in [0,2]$, (b) 2D map in spherical coordinates $(\psi_1, \psi_2)$ of $f(X,p)$ on the unit-sphere for $p = p_{k,fl}^*$.</figcaption></figure>

## 5. Conclusion

After having recalled the ki.s.s. issue and its solution for divergence of conservative and of non conservative systems as well, the 3 dimensional flutter ki.s.s. issue is investigated. First, using the usual algebraic criterion of flutter instability involving the discriminant of the characteristic polynomial, a sufficient condition involving the eigenvalues of the symmetric and the skew-symmetric parts of the operator is proposed that ensures a non optimal conditional flutter ki.s.s. characterized by the value $p_{s,k,fl}^*$ of the load parameter.

For getting necessary and sufficient conditions more advanced calculations are led. These conditions are brought back to an original minimization problem on the Grassmann manifold of the 2-planes $Gr_2(E)$. Two ways are used to tackle this issue. The first one uses differential calculation on this Grassmann manifold viewed as the base space of the principal fiber bundle of the Stiefel manifold $St_2(E)$ built by the families of two orthonormal vectors of $E$. It uses a more abstract point of view but leads to more compact calculations that may be generalized in higher dimensions. The second one uses the opportunity, in this three dimensional case, to parametrize the issue by the sphere because it is here a double cover of the involved Grassmann manifold $Gr_2(E)$. Calculations on the sphere are more usual and lead, after more tedious calculations, to the same results of a family of three explicit analytic nonlinear equations. The solutions of these equations lead to the optimal flutter ki.s.s. critical load $p_{k,fl}^*$ and to the corresponding kinematic constraint.

The mechanical consequences are then highlighted showing only a partial conditional flutter ki.s.s. understood in the sense of the safety of the structure. So, contrary to the divergence ki.s.s. – for the conservative and for the non conservative systems – the optimal flutter ki.s.s. condition does not ensure the flutter stability of the initial free system $\Sigma = \Sigma_{free}$ and the both cases $p_{k,fl}^* > p_{fl}^*$ or $p_{k,fl}^* < p_{fl}^*$ are possible. Finally, an illustration of these mechanical insights is done through a three degree of freedom Ziegler column with two distinct mass distributions.

## Acknowledgments

We would like to take the opportunity of this work to pay warm tribute to Professor Claude Vallée who very recently died. Since several decades, whatever the subject matter in theoretical mechanics (Lie group actions, dynamics, symmetry groups, rational mechanics, ...) the talks with Claude Vallée were always enlightening. We are personally sure that his original and fruitful contributions to mechanics will continue through his students and his friends and colleagues.

## References

[1] J. Lerbet, N. Challamel, F. Nicot, F. Darve, Variational formulation of divergence stability for constrained systems, Appl. Math. Model. (2015) http://dx.doi.org/10.1016/j.apm.2015.02.052.
[2] N. Challamel, F. Nicot, J. Lerbet, F. Darve, Stability of non-conservative elastic structures under additional kinematic constraints, Eng. Struct. 32 (2010) 3086-3092.

[3] J. Lerbet, M. Aldowaji, N. Challamel, F. Nicot, F. Prunier, F. Darve, P-positive definite matrices and stability of nonconservative systems, ZAMM Z. Angew. Math. Mech. 92 (5) (2012) 409-422.

[4] W. Thomson, P.G. Tait, Treatise on Natural Philosophy, in: Clarendon Press Series, 1879.

[5] T. Tarnai, *Paradoxical behaviour of vibrating systems challenging Rayleigh's theorem*, in: 21st International Congress of Theoretical and Applied Mechanics, Warsaw, 2004.

[6] J.M.T. Thompson, 'Paradoxical' mechanics under fluid flow, Nature 296 (5853) (1982) 135-137.

[7] A. Edelman, T.A. Arias, S. Smith, The geometry of algorithms with orthogonality constraints, SIAM J. Matrix Anal. Appl. 20 (2) (1998) 303-353.

[8] P.-A. Absil, R. Mahony, R. Sepulchre, Riemannian geometry of Grassmann manifolds with a view on algorithmic computation, Acta Appl. Math. 80 (2) (2004) 199-220.

[9] K.E. Gustafson, D.K.M. Rao, Numerical Range. The Field of Values of Linear Operators and Matrices, in: Universitext, Springer, 1997.

[10] W. Forst, D. Hoffmann, *Optimization—Theory and Practice*, Springer, 2010.