# Un modèle d'ordre réduit stable pour les problèmes vibro-acoustiques transitoires discrétisés par la méthode des éléments finis avec la formulation en $u - p - \phi$

L. Khoun<fnref n="1" />,<fnref n="2" />, C. Leblond<fnref n="2" />, M. Abbas<fnref n="3" />,<fnref n="4" />,
Y. Maday<fnref n="1" />, N. Aissiouene <fnref n="1" />, Y. Tampango <fnref n="3" />

**Résumé** — Dans ce travail, nous proposons une technique de construction de modèle d'ordre réduit qui préserve la stabilité du modèle complet dans les problèmes vibro-acoustiques transitoires formulés en $u - p - \phi$ (déplacement structurel, pression acoustique, potentiel de déplacement fluide) discrétisés avec la méthode des éléments finis. La technique proposée est également valable pour les problèmes extérieurs lorsque la condition de rayonnement de Sommerfeld par la méthode de BGT-1 [2] est utilisée. La stabilité du modèle d'ordre réduit est démontrée théoriquement et numériquement.

**Mots clés** — Modèle d'ordre réduit, préservation de la stabilité, vibro-acoustique transitoire.

## 1 Introduction

Malgré l'amélioration constante de la puissance de calcul des ordinateurs, le calcul et l'analyse de la réponse vibratoire de structures immergées demeurent des processus coûteux dans les industries navale et nucléaire, de par la taille et la complexité des structures en jeu, ainsi que leur couplage fort avec un fluide lourd. La discrétisation par éléments finis de tels problèmes vibro-acoustiques résulte en des systèmes matriciels à plusieurs millions de degrés de liberté, couplant les inconnues structure et fluide. Même pour les régimes harmoniques linéaires, le calcul de la réponse fréquentielle dans ce cadre est conséquent : les matrices à inverser à chaque fréquence sont complexes, non-hermitiennes et les systèmes résultants sont souvent mal conditionnés. Le temps de calcul total pour une réponse en bande fréquentielle fine peut dès lors devenir prohibitif en bureau d'études.

Les techniques de réduction de modèle sont maintenant capables de réduire drastiquement les temps de calcul sans perdre de précision. On peut trouver dans la littérature différentes méthodes de réduction, par exemple la Méthode des Bases Réduites (RBM) [3, 12], la POD [8], la PGD [4, 15] et la méthode de projection de Krylov [9, 10, 19].

Alors que la réduction du problème harmonique linéaire paramétrique [12] est mature pour certaines applications industrielles, tels les problèmes de discrétion acoustique [13], d'autres applications industrielles requièrent le postraitement d'historique temporel des contraintes, et donc nécessitent la réduction du problème vibro-acoustique transitoire paramétrique.

Pour les problèmes transitoires, il est nécessaire que le modèle réduit ait la même propriété de stabilité que le modèle complet discrétisé. Dans les problèmes vibro-acoustiques transitoires intérieurs (i.e le domaine fluide est fini) formulés en $u - p$ (déplacement structurel - pression acoustique) et $u - \psi$ (déplacement structurel - potentiel de vitesse de déplacement fluide) discrétisés par la méthode des éléments finis, une technique de construction de modèle d'ordre réduit qui préserve la stabilité a été proposée dans [18]. Cette technique a été ensuite déployée dans [16] pour les problèmes extérieurs (i.e le domaine fluide est infini) qui utilisent la méthode des méthodes des éléments finis [1] pour modéliser la condition de rayonnement de Sommerfeld.

Dans ce travail, nous adaptons la technique proposée dans [18] pour construire des modèles réduits stables pour les problèmes vibro-acoustiques transitoires formulés en $u - p - \phi$. Pour les problèmes extérieurs, la technique proposée restera valable si on modélise la condition de rayonnement de Sommerfeld

<footnote n="1">¹ Laboratoire J.-L. Lions, Sorbonne Université, ladya.khoun@etu.upmc.fr, {yvon.maday, nora.aissiouene}@upmc.fr</footnote>

<footnote n="2">² Naval Group Research, {cedric.leblond, ladya.khoun} @naval-group.com</footnote>

<footnote n="3">³ EDF R & D {mickael.abbas, yannick.tampango}@edf.fr</footnote>

<footnote n="4">⁴ Institut des Sciences de la Mécanique et Applications Industrielles, UMR 9219 EDF-CNRS-CEA-ENSTA, Université Paris Saclay</footnote>
## 2 Modèle complet et la stabilité

Le modèle complet se réfère au modèle discret obtenu par la discrétisation des variables d'espaces par la méthode éléments finis. Le modèle complet est gouverné par un système d'EDO d'ordre 2 :

$$ egin{cases} M\ddot{X}(t) + C\dot{X}(t) + KX(t) = f(t)F \\ X(0) = X_0, \dot{X}(0) = \dot{X}_0 \end{cases} \quad (1) $$

où $M, C$ et $K \in \mathbb{R}^{n \times n}$ sont respectivement la matrice de masse, d'amortissement et de rigidité, $F \in \mathbb{R}^n$ est le vecteur second membre du problème et $f$ est une fonction dépendant du temps. Le modèle (1) est dit **stable** si les racines du polynôme $P(s) = \det(s^2M + sC + K)$ sont toutes de partie réelle négative.

Des différentes formulations en vibro-acoustique peut être trouvées dans la littérature [5, 6, 7, 14, 17]. Nous nous intéressons dans ce travail à la formulation en $u - p - \phi$ qui s'écrit :

$$ \underbrace{\begin{pmatrix} M_s & 0 & -\rho_f K_c \\ 0 & 0 & \frac{1}{c_f^2} M_f \\ -\rho_f K_c^T & \frac{1}{c_f^2} M_f & -\frac{\rho_f}{c_f^2} K_f \end{pmatrix}}_{M_{up\phi}} \begin{pmatrix} \ddot{U}(t) \\ \ddot{P}(t) \\ \ddot{\Phi}(t) \end{pmatrix} + \underbrace{\begin{pmatrix} C_s & 0 & 0 \\ 0 & 0 & 0 \\ 0 & C_f & 0 \end{pmatrix}}_{C_{up\phi}} \begin{pmatrix} \dot{U}(t) \\ \dot{P}(t) \\ \dot{\Phi}(t) \end{pmatrix} + \underbrace{\begin{pmatrix} K_s & 0 & 0 \\ 0 & \frac{1}{\rho_f c_f^2} M_f & 0 \\ 0 & 0 & 0 \end{pmatrix}}_{K_{up\phi}} \begin{pmatrix} U(t) \\ P(t) \\ \Phi(t) \end{pmatrix} = f(t) F_{up\phi} \quad (2) $$

où $\rho_f$ est la masse volumique du fluide, $c_f$ est la célérité du son dans le fluide, $M_s, C_s$ et $K_s$ sont respectivement la matrice de masse, d'amortissement et de rigidité de la structure, $M_f, C_f$ et $K_f$ sont respectivement la matrice de masse, d'amortissement et de rigidité du fluide et $K_c$ la matrice de couplage. La matrice $M_s, M_f, K_s, K_f$ sont symétriques définies positives. Nous nous limitons dans le cas où la matrice d'amortissement $C_s$ et $C_f$ sont symétriques semi-définies positives.

## 3 Modèle d'ordre réduit stable basé sur la projection de Galerkin

### 3.1 Construction du modèle réduit

Soit $V \in \mathbb{C}^{n \times N}$, nous nous intéresserons au modèle réduit obtenu par la projection de Galerkin sur la base $V$ qui s'écrit :

$$ egin{cases} M_r \ddot{X}_r(t) + C_r \dot{X}_r(t) + K_r X_r(t) = f(t) F_r \\ X_r(0) = V^* X_0, \dot{X}_r(0) = V^* \dot{X}_0 \end{cases} \quad (3) $$

avec $M_r = V^*MV, C_r = V^*CV$ et $K_r = V^*KV \in \mathbb{C}^{N \times N}$ et $F_r = V^*F \in \mathbb{C}^N$. L'approximation de la solution du modèle complet (1) par le modèle réduit (3) est : $X^{MOR}(t) = V X_r(t)$. En général, le modèle réduit (3) et le modèle complet (1) n'admettent pas la même propriété de stabilité.

Nous proposons ici de construire les modèles d'ordre réduit par la projection de Galerkin sur la base $V$ de la forme :

$$ V = \begin{pmatrix} V_u & 0 & 0 \\ 0 & V_p & 0 \\ 0 & 0 & V_p \end{pmatrix} \in \mathbb{C}^{n \times (N_u + 2N_p)} \quad (4) $$

afin de préserver la stabilité du modèle complet (2). En effet, nous pouvons montrer, par le calcul du polynôme $P(s)$ défini dans la section précédente, que le modèle réduit ainsi construit admet la même propriété de stabilité que le modèle réduit obtenu par la projection de Galerkin du modèle complet en $u - p$ sur la base $V_{up} = \begin{pmatrix} V_u & 0 \\ 0 & V_p \end{pmatrix} \in \mathbb{C}^{n \times (N_u + N_p)}$ [11]. Dans [18], il a été démontré que dans la formulation en $u - p$, le modèle d'ordre réduit obtenu par la projection de Galerkin sur la base $V_{up}$ est stable. Nous
pouvons donc conclure que dans la formulation en $u - p - \phi$, le modèle d'ordre réduit obtenu par la projection de Garlerkin sur la base de la forme (4) est stable.

La traduction de la condition de rayonnement de Sommerfeld par la méthode BGT-1 [2] pour traiter le problème du domaine fluide infini équivalant à ajouter une matrice d'amortissement symétrique semi-définie positive $C_f^{BGT}$ dans le fluide et modifier la matrice de rigidité du fluide $K_f$ par $K_f^{new} = K_f + K_f^{BGT}$ où $K_f^{BGT}$ est une matrice symétrique semi-définie positive. La technique proposée reste par conséquent valable dans ce cas.

## 3.2 Construction de la base réduite

La qualité du modèle d'ordre réduit obtenu par la projection de Garlerkin présenté dans la section 3.1 dépend fortement du choix de la base réduite $V$. Différentes méthodes de construction de la base réduite $V$ peuvent être retrouvées dans la littérature [8, 9, 10, 19]. Nous proposons ici d'utiliser la base réduite obtenue par l'algorithme glouton appliqué sur le problème complet (2) en régime harmonique en considérant la fréquence comme le paramètre du problème. Le problème complet (2) en régime harmonique s'écrit :

$$ \underbrace{[-\omega^2 M_{up\phi} + i\omega C_{up\phi} + K_{up\phi}]}_{A_{up\phi}(\omega)} \tilde{X}_{up\phi}(\omega) = F_{up\phi} \quad (5) $$

où $\omega$ représente la pulsation. La construction classique de la base réduite par l'algorithme glouton est présentée dans l'algorithme suivant :

**Algorithm 1 :** Construction de la base réduite avec l'algorithme glouton

**Entrée** : $N_{max}$ (nombre de mode), $N_{tir}$ (nombre d'échantillonnages à tirer aléatoirement à chaque itération), $\omega_{min}$ (la pulsation minimale d'intérêt), $\omega_{max}$ (la pulsation maximale d'intérêt), les matrices $M_{up\phi}, C_{up\phi}, K_{up\phi}$ et le second membre $F_{up\phi}$.

**Sortie** : La base réduite $V_{N_{max}}$

1. On choisit $\tilde{X}_{up\phi}(\omega_{min})$ pour le premier mode :

2. Enrichir la base :

$$ V_1 = \left\{ \frac{\tilde{X}_{up\phi}(\omega_{min})}{\|\tilde{X}_{up\phi}(\omega_{min})\|} \right\} $$

3. **while** ($k \ge 2$ and $k \le N_{max}$) **do**

4. Générer un échantillon aléatoire $S_{N_{tir}}$ de $N_{tir}$ valeurs dans $[\omega_{min}, \omega_{max}]$

5. Résoudre pour chaque valeur $\omega \in S_{N_{tir}}$ le problème réduit de taille $k-1$ :

$$ A_{up\phi}^r(\omega) \tilde{X}_{up\phi}^r(\omega) = F_{up\phi}^r $$

$$ \text{où } A_{up\phi}^r(\omega) = V_{k-1}^* A_{up\phi}(\omega) V_{k-1} \text{ et } F_{up\phi}^r = V_{k-1}^* F_{up\phi}. $$

6. Identifier la valeur $\omega^*$ qui maximise la norme de résidus :

$$ \omega^* = \underset{\omega \in S_{N_{tir}}}{\operatorname{argmax}} \|F_{up\phi} - A_{up\phi}(\omega)V_{k-1}\tilde{X}_{up\phi}^r(\omega)\| $$

7. Orthogonaliser le nouveau mode $\tilde{X}_{up\phi}(\omega^*)$ avec la base $V_{k-1}$

8. Enrichir la base :

$$ V_k = [V_{k-1}, \text{ortho}(\tilde{X}_{up\phi}(\omega^*), V_{k-1})] $$

9. $k = k + 1$

10. **end while**

Afin de s'assurer la stabilité du modèle réduit, nous proposons de transformer la base au fur et à mesure à la forme (4) dans la boucle de l'algorithme glouton. Le nouvel algorithme glouton est présenté dans l'algorithme 2.
**Algorithm 2** : Construction de la base réduite qui assure la stabilité du modèle d'ordre réduit avec l'algorithme glouton

**Entrée** : $N_{max}$ (nombre de mode), $N_{tir}$ (nombre d'échantillonnages à tirer aléatoirement à chaque itération), $\omega_{min}$ (la pulsation minimale d'intérêt), $\omega_{max}$ (la pulsation maximale d'intérêt), les matrices $M_{up\phi}, C_{up\phi}, K_{up\phi}$ et le second membre $F_{up\phi}$.

**Sortie** : La base réduite $V_{N_{max}}$

1: On choisit $\tilde{X}_{up\phi}(\omega_{min}) = (U(\omega_{min}), P(\omega_{min}), \Phi(\omega_{min}))^T$ pour le premier mode

2: Enrichir la base :

$$V_1 = \begin{pmatrix} \frac{U(\omega_{min})}{\|U(\omega_{min})\|} & 0 & 0 \\ 0 & \frac{P(\omega_{min})}{\|P(\omega_{min})\|} & 0 \\ 0 & 0 & \frac{P(\omega_{min})}{\|P(\omega_{min})\|} \end{pmatrix}$$

3: **while** ($k \ge 2$ and $k \le N_{max}$) **do**

4: Générer un échantillon aléatoire $S_{N_{tir}}$ de $N_{tir}$ valeurs dans $[\omega_{min}, \omega_{max}]$

5: Résoudre pour chaque valeur $\omega \in S_{N_{tir}}$ le problème réduit de taille $3(k-1)$ :

$$A_{up\phi}^r(\omega) \tilde{X}_{up\phi}^r(\omega) = F_{up\phi}^r$$

$$\text{où } A_{up\phi}^r(\omega) = V_{k-1}^* A_{up\phi}(\omega) V_{k-1} \text{ et } F_{up\phi}^r = V_{k-1}^* F_{up\phi}.$$

6: Identifier la valeur $\omega^*$ qui maximise la norme de résidus :

$$\omega^* = \operatorname*{argmax}_{\omega \in S_{N_{tir}}} \| F_{up\phi} - A_{up\phi}(\omega) V_{k-1} \tilde{X}_{up\phi}^r(\omega) \|$$

7: Orthogonaliser les vecteurs $(U(\omega^*), 0, 0)^T, (0, P(\omega^*), 0)^T, (0, 0, P(\omega^*))^T$ avec la base $V_{k-1}$

8: Enrichir la base :

$$V_k = [V_{k-1}, \text{ortho}((U(\omega^*), 0, 0)^T, V_{k-1}), \text{ortho}((0, P(\omega^*), 0)^T, V_{k-1}), \text{ortho}((0, 0, P(\omega^*))^T, V_{k-1})]$$

9: $k = k + 1$

10: **end while**

# 4 Validation numérique

## 4.1 Cas d'étude

La géométrie du cas d'étude numérique est présentée par la figure 1. Il s'agit d'une coque circulaire cylindrique de longueur finie, et d'une plaque encastrée à ses extrémités, immergées dans un fluide infini. La coque cylindrique est de longueur $0.5 \, m$, de rayon $0.1 \, m$ et de l'épaisseur $0.01 \, m$. La plaque est carrée de côté $0.5 \, m$ et d'épaisseur $0.01 \, m$. Pour modéliser la condition de rayonnement de Sommerfeld, le domaine fluide est tronqué en une demie-sphère de rayon $1 \, m$ et la condition de BGT-1 est appliquée sur le bord sphérique du domaine fluide. Le fluide est de masse volumique $\rho_f = 1000 \, kg \cdot m^{-3}$. La célérité du son dans le fluide est $c_f = 1500 \, m \cdot s^{-1}$. Les matériaux constituant la structure sont tous de comportement élastique linéaire isotrope avec le module Young $E = 2.1 \cdot 10^{11} \, Pa$, le coefficient de Poisson $\nu = 0.3$ et la masse volumique $\rho_s = 7850 \, kg \cdot m^{-3}$. La sollicitation est introduite par un effort nodal $f(t) \in_x$ sur le nœud positionné au point $(0.1, 0, 0.5)$ (voir la figure 1).

Pour la discrétisation du problème par la méthode des éléments finis, on utilise les éléments plaques DKT [5] pour la structure et les éléments tétraèdres linéaires pour le fluide. Le modèle d'amortissement de Rayleigh est utilisé pour modéliser l'amortissement dans la structure ($C_s = \alpha_s K_s$ avec $\alpha_s = 10^{-5}$). Le modèle complet dans notre étude est un système d'EDO d'ordre 2 de taille $n = 32298$ :

$$\left\{ \begin{array}{l} M_{up\phi} \ddot{X}_{up\phi}(t) + C_{up\phi} \dot{X}_{up\phi}(t) + K_{up\phi} X_{up\phi}(t) = f(t) F_{up\phi} \\ X_{up\phi}(0) = 0, \dot{X}_{up\phi}(0) = 0 \end{array} \right. \quad (6)$$

où $F_{up\phi} \in \mathbb{R}^n$ est le vecteur second membre, correspondant à la force nodal unité sur le nœud positionné au
<figure><img src="image_1.png" /><figcaption>FIGURE 1 – Illustration du maillage du cas d'étude.</figcaption></figure>

point (0.1, 0, 0.5). Pour que les matrices soient mieux conditionnées, l'ensemble des données physiques du problème est adimensionné à partir des échelles caractéristiques fournies par le tableau 1.

<table><caption>TABLE 1 – Échelles caractéristiques du problème adimensionné</caption><thead><tr><td>Variable</td><td>Échelles caractéristiques</td></tr></thead><tbody><tr><td>Longueur, déplacement</td><td>L = 1m</td></tr><tr><td>Temps</td><td>L/c<sub>f</sub></td></tr><tr><td>Fréquence</td><td>c<sub>f</sub>/L</td></tr><tr><td>Masse volumique</td><td>ρ<sub>f</sub></td></tr><tr><td>Pression, contrainte</td><td>ρ<sub>f</sub>c<sub>f</sub><sup>2</sup></td></tr><tr><td>Potentiel</td><td>L<sup>2</sup></td></tr><tr><td>Impédance</td><td>ρ<sub>f</sub>c<sub>f</sub></td></tr></tbody></table>

## 4.2 Illustration de la stabilité du modèle réduit

Pour vérifier numériquement la stabilité du modèle d'ordre réduit, nous pouvons calculer leurs pôles pour apprécier s'ils sont tous de parties réelles négative. La partie réelle non-négative la plus grande des pôles du modèle d'ordre réduit obtenu par la projection de Galerkin sur la base issue de l'algorithme 1 et 2 est présentée par la figure 2.

<figure><img src="image_2.png" /><figcaption>FIGURE 2 – La partie réelle non-négative la plus grande des pôles du modèle d'ordre réduit obtenu par la projection de Galerkin sur la base issue de l'algorithme 1 (à gauche) et l'algorithme 2 (à droite) en fonction de la dimension du modèle d'ordre réduit.</figcaption></figure>
D'après la figure 2, nous observons que le modèle d'ordre réduit obtenu par la projection de Galerkin sur la base issue de l'algorithme 1 n'est pas toujours stable. Il y a des instabilités dès N = 8 pour notre cas test. Pour les modèles d'ordre réduit construits par la base issue de l'algorithme 2, nous voyons clairement qu'ils sont toujours stables. De plus, nous voyons dans la figure 3 que la solution du modèle d'ordre réduit construit par la projection de Galerkin sur la base issue de l'algorithme 2 fournit une très bonne approximation de la solution du modèle complet.

<figure><img src="image_3.png" /><figcaption>FIGURE 3 – La pression moyenne sur la surface du cylindre obtenue par le modèle complet et le modèle d'ordre réduit construit par la projection de Galerkin sur la base issue de l'algorithme 2 pour $N_{max} = 15$. À gauche, cas $f(t) = e^{-(a/t+bt)}\sin(\frac{2\pi t}{T})$ avec $a = 0.1, b = 2.0$ et $T = 0.2$. À droite, cas $f(t) = e^{-(a/t+bt)}\sin(\frac{2\pi t}{T})$ avec $a = 0.1, b = 2.5$ et $T = 0.5$.</figcaption></figure>

# 5 Conclusion

Dans ce papier, une technique de construction de modèle réduit qui préserve la stabilité du modèle complet discret issue de la discrétisation éléments finis de la formulation en $u - p - \phi$ est proposée et validée numériquement. La technique de construction de la base réduite qui garantit de préserver la stabilité du modèle complet s'appuie sur un algorithme glouton. L'étude numérique atteste que le modèle réduit construit par la projection de Galerkin sur cette base est à la fois stable et précis.

# Références

[1] R. J. Astley, J. P. Coyette, L. Cremers, *Three-dimensional wave-envelope elements of variable order for acoustic radiation and scattering. Part II. Formulation in the time domain*, J. Acoust. Soc. Am. 103 (1) (1998) 64-72.

[2] A. Bayliss, M. Gunzburger, E. Turkel, *Boundary conditions for the numerical solution of elliptic equations in exterior regions*, SIAM Journal on Applied Mathematics, vol. 42, 430-451, 1982.

[3] Y. Chen, J. S. Hesthaven, Y. Maday, J. Rodriguez, *Certified reduced basis methods and output bounds for the harmonic Maxwell's equations*, SIAM J. Sci. Comput., 32(2), 970-996.

[4] A. Dumon, C. Allery, A. Ammar, *Proper Generalized Decomposition method for incompressible Navier-Stokes equations with a spectral discretization*, Applied Mathematics and Computation 219 (2013) 8145-8162.

[5] Électricité de France, *Finite element code_aster, Structures and Thermomechanics Analysis for Studies and Research*, Open source on www.code-aster.org, 1989-2018.

[6] G. C. Everstine, *Finite element formulatons of structural acoustics problems*, Computers & Structure Vol. 65, No.3, pp. 307-321, 1997.

[7] G. C. Everstine, *A symmetric potentiel formulation for fluide-structure interaction*, Journal of Sound and Vibration (1981) 79(1), 157-166.

[8] L. Fick, Y. Maday, A. T. Patera, T. Taddei, *A Reduced Basis Technique for Long-Time Unsteady Turbulent Flows*, Journal of Computational Physics, October 2017.

[9] E. J. Grimme, *Krylov Projection Methods for Model Reduction*, PhD thesis, Univ. of Illinois, Urbana-Champaign, (1997).
[10] S. Gugercin, A. C. Antoulas, *Rational Krylov Methods for Optimal H2 Model Reduction* (2006).

[11] L. Khoun, *Rapport technique interne Naval Group/EDF/Laboratoire Jacques-Louis Lions* (2018).

[12] C. Leblond, J.-F. Sigrist, *A reduced basis approach for the parametric low frequency response of submerged viscoelastic structures*, Finite Elements in Analysis and Design 119 (2016) 15-29.

[13] C. Leblond, M. Chevreuil, C. Allery, C. Beghein, *A goal-oriented model reduction technique for parametric fluid-structure problems*, ECCMG, Glasgow, UK, June 2018.

[14] H. Morand, R. Ohayon, *Fluid Structure Interaction: Applied Numerical Methods*, Wiley, 1995.

[15] A. Nouy, *A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations*, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2010, 199 (23-24), pp.1603-1626.

[16] S.van Ophem, O. Atak, E. Deckers, W. Desmet, *Stable model order reduction for time-domain exterior vibro-acoustic finite element simulations*, CMAME 325 (2017), 240-264.

[17] J.-F. Sigrist, *Fluid-Structure Interaction: An Introduction to Finite Element Coupling*, Wiley, 2015.

[18] A. van de Walle, F. Naets, E. Deckers, W. Desmet, *Stability-preserving model-order reduction for time-domain simulation of vibro-acoustic FE models*, Internat. J. Numer. Methods Engrg. 109 (6) (2017) 889-912.

[19] S. Wyatt, *Issues in Interpolatory Model Reduction: Inexact Solves, Second-order Systems and DAEs*, PhD thesis, Univ. Virginia, (2012).