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

HAL
open science

# Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime

Lucie Gallerand, Mathias Legrand, Raymond Panneton, Philippe Leclaire, Thomas Dupont

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

## To cite this version:

Lucie Gallerand, Mathias Legrand, Raymond Panneton, Philippe Leclaire, Thomas Dupont. Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime. 2024. [https://hal.science/hal-04528668](https://hal.science/hal-04528668)

HAL Id: hal-04528668

[https://hal.science/hal-04528668](https://hal.science/hal-04528668)

Preprint submitted on 30 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.
# Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime

Lucie Gallerand, Mathias Legrand, Raymond Panneton, Philippe Leclaire,
Thomas Dupont

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

## To cite this version:

Lucie Gallerand, Mathias Legrand, Raymond Panneton, Philippe Leclaire, Thomas Dupont. Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime. Journal of Sound and Vibration, 2025, pp.119018. 10.1016/j.jsv.2025.119018. hal-04528668v2

HAL Id: hal-04528668

https://hal.science/hal-04528668v2

Submitted on 14 Mar 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime

Lucie Gallerand<sup>a</sup>, Mathias Legrand<sup>b</sup>, Raymond Panneton<sup>c</sup>, Philippe Leclaire<sup>d</sup>, Thomas Dupont<sup>b,a</sup>

<sup>a</sup>Department of Mechanical Engineering, École de technologie supérieure, Montréal, Canada

<sup>b</sup>Department of Mechanical Engineering, McGill University, Montréal, Canada

<sup>c</sup>CRASH-UdeS, Department of Mechanical Engineering, Université de Sherbrooke, Sherbrooke, Canada

<sup>d</sup>DRIVE EA1859, Université de Bourgogne Franche-Comté, ISAT, Nevers, France

Microperforated plates (MPP) are commonly used in the field of acoustics for sound absorption purposes. Recent research in the area of structural dynamics revealed that they also provide additional viscous damping in the low-frequency range because of fluid-solid interactions in the microperforation boundary layers, mainly through viscous friction mechanisms. It is now established that the vibratory behavior of these systems can be modelled using a homogenization procedure that produces a pair of coupled partial differential equations (PDEs) collectively governing the dynamics of both a structural plate and a fictitious virtual fluid plate. It has been observed that the added damping achieves its maximum at a characteristic frequency governed in particular by the size of the perforations. It is also known that these systems, which can be implemented in hostile environments such as aircraft turbines, involve nonlinear mechanisms for large mechanical and/or acoustic excitations. Two causes of nonlinearities can be distinguished: (i) one associated with a high amplitude harmonic speed of the fluid in the perforations, and (ii) one associated with large harmonic structural deformations. The two causes can combine under strong excitation. The present paper focuses exclusively on the first cause of nonlinearity. As the fluid velocity increases inside the perforations, resistive and inertial phenomena arise within the perforations, influencing the structural response. These effects can be partly captured analytically by the Forchheimer correction, materialized by an additional quadratic nonlinear damping term in the governing PDE governing the dynamics of the fictitious fluid plate. The governing equations are solved numerically, and sensitivity analyses are carried out on the added damping to the excitation level. The proposed model is validated by experiments conducted on an equivalent cantilevered MPP. Analytical and experimental results show that the added viscous damping depends on the relative fluid-solid velocity. The added damping effect can, depending on the perforation diameter, reach a maximum for a critical value of the relative fluid-solid velocity, with all other independent parameters fixed. In the nonlinear framework, the added damping is a function of space and depends not only on the perforation diameter as in the linear framework, but also on the relative fluid-solid velocity, and is defined as a function of space.

Keywords: Microperforated plates — Added nonlinear damping — Forchheimer’s law

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Model</td><td>4</td></tr><tr><td>2.1</td><td>Dynamics of a finite-size microperforated plate saturated by a light fluid</td><td>4</td></tr><tr><td>2.2</td><td>Added damping in linear regime</td><td>5</td></tr><tr><td>2.3</td><td>Added damping in nonlinear regime</td><td>5</td></tr><tr><td>3</td><td>Resolution procedure</td><td>6</td></tr><tr><td>3.1</td><td>Harmonic method balance</td><td>6</td></tr><tr><td>3.2</td><td>Convergence analysis</td><td>8</td></tr><tr><td>4</td><td>Sensitivity of microperforated plates vibratory response to nonlinear damping</td><td>9</td></tr><tr><td>4.1</td><td>Frequency response</td><td>9</td></tr><tr><td>4.2</td><td>Normalized resistance and maximum of added damping in nonlinear regime</td><td>10</td></tr><tr><td>4.3</td><td>From a characteristic frequency to a critical relative velocity</td><td>11</td></tr><tr><td>4.4</td><td>Space distribution of the added damping</td><td>12</td></tr><tr><td>5</td><td>Experimental validation</td><td>14</td></tr><tr><td>5.1</td><td>Experimental set-up</td><td>14</td></tr><tr><td>5.2</td><td>Relative displacement</td><td>14</td></tr><tr><td>5.3</td><td>Experimental and analytical results</td><td>15</td></tr><tr><td>6</td><td>Conclusions</td><td>16</td></tr><tr><td></td><td>References</td><td>17</td></tr></table>

Email addresses: lucie.gallerand.1@ens.etsmtl.ca (Lucie Gallerand), mathias.legrand@mcgill.ca (Mathias Legrand), raymond.panneton@USherbrooke.ca (Raymond Panneton), philippe.leclaire@u-bourgogne.fr (Philippe Leclaire), thomas.dupont@etsmtl.ca (Thomas Dupont)
# Added nonlinear damping of homogenized fluid-saturated microperforated plates in Forchheimer flow regime

Lucie Gallerand<sup>a,*</sup>, Mathias Legrand<sup>b</sup>, Raymond Panneton<sup>c</sup>, Philippe Leclaire<sup>d</sup>, Thomas Dupont<sup>e,a</sup>

Microperforated plates (MPP) are commonly used in the field of acoustics for sound absorption purposes. Recent research in the area of structural dynamics revealed that they also provide additional viscous damping in the low-frequency range because of fluid-solid interactions in the microperforation boundary layers, mainly through viscous friction mechanisms. It is now established that the vibratory behavior of these systems can be modelled using a homogenization procedure that produces a pair of coupled partial differential equations (PDEs) collectively governing the dynamics of both a structural plate and a fictitious virtual fluid plate. It has been observed that the added damping achieves its maximum at a characteristic frequency governed in particular by the size of the perforations. It is also known that these systems, which can be implemented in hostile environments such as aircraft turbines, involve nonlinear mechanisms for large mechanical and/or acoustic excitations. Two causes of nonlinearities can be distinguished: (i) one associated with a high amplitude harmonic speed of the fluid in the perforations, and (ii) one associated with large harmonic structural deformations. The two causes can combine under strong excitation. The present paper focuses exclusively on the first cause of nonlinearity. As the fluid velocity increases inside the perforations, resistive and inertial phenomena arise within the perforations, influencing the structural response. These effects can be partly captured analytically by the Forchheimer correction, materialized by an additional quadratic nonlinear damping term in the PDE governing the dynamics of the fictitious fluid plate. The system of equations is solved numerically, and sensitivity analyses are carried out on the added damping to the excitation level. The proposed model is validated by experiments conducted on an equivalent cantilevered MPP. Analytical and experimental results show that the added viscous damping depends on the relative fluid-solid velocity. The added damping effect can, depending on the perforation diameter, reach a maximum for a critical value of the relative fluid-solid velocity, with all other independent parameters fixed. In the nonlinear framework, the added damping is a function of space and depends not only on the perforation diameter as in the linear framework, but also on the relative fluid-solid velocity, and is defined as a function of space.

**Keywords:** Microperforated plates — Added nonlinear damping — Forchheimer's law

## Contents

1 Introduction 3

2 Model 4

2.1 Dynamics of a finite-size microperforated plate saturated by a light fluid 4

2.2 Added damping in linear regime 5

2.3 Added damping in nonlinear regime 6

3 Solution method 6

3.1 Harmonic balance method 6

3.2 Convergence analysis 8

4 Sensitivity of microperforated plates vibratory response to nonlinear damping 9

4.1 Frequency response 9

4.2 Maximum of added damping in nonlinear regime 10

4.3 Space distribution of the added damping 12

<footnote><sup>a</sup>Department of Mechanical Engineering, École de technologie supérieure, 1100 Rue Notre Dame O, Montréal, QC H3C1K3, Canada</footnote>

<footnote><sup>b</sup>Department of Mechanical Engineering, McGill University, 845 Rue Sherbrooke O, Montréal, QC H3A0G4, Canada</footnote>

<footnote><sup>c</sup>CRASH-UdeS, Department of Mechanical Engineering, Université de Sherbrooke, 2500 Bd de l'Université, Sherbrooke, QC J1N3C6, Canada</footnote>

<footnote><sup>d</sup>DRIVE EA1859, Université de Bourgogne Franche-Comté, ISAT, 49 Rue Mademoiselle Bourgeois, 58027 Nevers, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: lucie.gallerand.1@ens.etsmtl.ca (Lucie Gallerand), mathias.legrand@mcgill.ca (Mathias Legrand), raymond.panneton@USherbrooke.ca (Raymond Panneton), philippe.leclaire@u-bourgogne.fr (Philippe Leclaire), thomas.dupont@etsmtl.ca (Thomas Dupont)</footnote>
# Nomenclature

Primary quantities introduced in the paper are recalled below with emphasis on the main explicit dependencies to other quantities. They are no repeated in the text. When deemed appropriate, secondary quantities with appropriate dependencies are directly introduced in the text.

## Geometric constants

<table><thead><tr><th>Symbol</th><th>Description</th><th>Explicit dependencies</th><th>Unit</th></tr></thead><tbody><tr><td>d</td><td>diameter of perforation</td><td>–</td><td>m</td></tr><tr><td>h</td><td>plate thickness</td><td>–</td><td>m</td></tr><tr><td>L<sub>x</sub></td><td>plate length in x direction</td><td>–</td><td>m</td></tr><tr><td>L<sub>y</sub></td><td>plate length in y direction</td><td>–</td><td>m</td></tr></tbody></table>

## Macroscopic parameters

<table><thead><tr><th>Symbol</th><th>Description</th><th>Explicit dependencies</th><th>Unit</th></tr></thead><tbody><tr><td>α<sub>∞</sub></td><td>effective tortuosity</td><td>φ, d</td><td>–</td></tr><tr><td>φ</td><td>perforation ratio</td><td>d</td><td>–</td></tr><tr><td>σ</td><td>nonlinear fluid-flow resistivity</td><td>φ, d, x, y, t</td><td>N s m<sup>-4</sup></td></tr><tr><td>σ<sub>0</sub></td><td>static fluid-flow resistivity</td><td>φ, d</td><td>N s m<sup>-4</sup></td></tr><tr><td>ε</td><td>Forchheimer's parameter</td><td>φ, d, h</td><td>s m<sup>-1</sup></td></tr><tr><td>C<sub>d</sub></td><td>discharge coefficient</td><td>φ, d, h</td><td>–</td></tr></tbody></table>

## Mechanical parameters

<table><thead><tr><th>Symbol</th><th>Description</th><th>Explicit dependencies</th><th>Unit</th></tr></thead><tbody><tr><td>α</td><td>Biot's elastic coupling coefficient</td><td>–</td><td>–</td></tr><tr><td>η</td><td>loss factor</td><td>–</td><td>–</td></tr><tr><td>η<sub>i</sub></td><td>nonlinear loss factor of mode i</td><td>φ, d, ω, ε, ḋ</td><td>–</td></tr><tr><td>η<sub>i</sub><sup>L</sup></td><td>linear loss factor of mode i</td><td>φ, d, ω</td><td>–</td></tr><tr><td>μ<sub>f</sub></td><td>fluid dynamic viscosity</td><td>–</td><td>kg m<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>ν</td><td>Poisson's coefficient</td><td>–</td><td>–</td></tr><tr><td>ρ</td><td>plate effective density</td><td>φ</td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>f</sub></td><td>fluid density</td><td>–</td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>s</sub></td><td>structure density</td><td>–</td><td>kg m<sup>-3</sup></td></tr><tr><td>D</td><td>bending stiffness</td><td>φ</td><td>GPa m<sup>3</sup></td></tr><tr><td>E</td><td>Young's modulus</td><td>–</td><td>GPa</td></tr><tr><td>K<sub>f</sub></td><td>fluid bulk modulus</td><td>–</td><td>kPa</td></tr><tr><td>M<sub>f</sub></td><td>Biot's elastic coefficient</td><td>–</td><td>N m<sup>-2</sup></td></tr></tbody></table>

## Other parameters

<table><thead><tr><th>Symbol</th><th>Description</th><th>Explicit dependencies</th><th>Unit</th></tr></thead><tbody><tr><td>(x, y)</td><td>coordinate system attached to plate</td><td>–</td><td>m</td></tr><tr><td>(x<sub>e</sub>, y<sub>e</sub>)</td><td>coordinate of excitation point</td><td>–</td><td>m</td></tr><tr><td>(x<sub>g</sub>, y<sub>g</sub>)</td><td>coordinate of observation point</td><td>–</td><td>m</td></tr><tr><td>x</td><td>vector of Fourier coefficients</td><td>ω</td><td>m</td></tr><tr><td>ω</td><td>driving angular frequency</td><td>–</td><td>rad s<sup>-1</sup></td></tr><tr><td>ω<sub>i</sub></td><td>MPP angular frequency of mode i</td><td>–</td><td>rad s<sup>-1</sup></td></tr><tr><td>Ψ<sub>i</sub></td><td>modal shape of mode i</td><td>x, y</td><td>–</td></tr><tr><td>f̃<sub>ext</sub></td><td>forcing force in frequency domain</td><td>ω</td><td>N</td></tr><tr><td>f̃<sub>NL</sub></td><td>nonlinear force in frequency domain</td><td>x, ω</td><td>N</td></tr><tr><td>f<sub>c</sub></td><td>characteristic frequency</td><td>φ, d</td><td>Hz</td></tr><tr><td>F<sub>ext</sub></td><td>forcing force amplitude</td><td>–</td><td>N</td></tr><tr><td>f<sub>ext</sub></td><td>forcing force in time domain</td><td>x, y, t, ω</td><td>N</td></tr><tr><td>f<sub>NL</sub></td><td>nonlinear force in time domain</td><td>x, y, t</td><td>N</td></tr><tr><td>N</td><td>number of degrees-of-freedom</td><td>–</td><td>–</td></tr><tr><td>N<sub>h</sub></td><td>number of harmonics</td><td>–</td><td>–</td></tr><tr><td>R<sub>L</sub></td><td>normalized linear resistance</td><td>φ, d</td><td>–</td></tr><tr><td>R<sub>NL</sub></td><td>normalized nonlinear resistance</td><td>φ, d, x, y, w</td><td>–</td></tr><tr><td>t</td><td>time</td><td>–</td><td>s</td></tr><tr><td>V</td><td>amplitude of ḋ</td><td>x, y, φ, d</td><td>m s<sup>-1</sup></td></tr><tr><td>V<sub>cr</sub></td><td>critical relative velocity</td><td>φ, d, ω</td><td>m s<sup>-1</sup></td></tr><tr><td>w</td><td>fluid-solid relative displacement</td><td>x, y, t</td><td>m</td></tr><tr><td>w<sub>s</sub></td><td>solid displacement</td><td>x, y, t</td><td>m</td></tr></tbody></table>

# 1. Introduction

Microperforated plates (MPP) are commonly used as lightweight acoustic-absorbing materials, serving as a viable substitute for porous materials and conventional acoustic resonators. MPP have versatile applications in acoustic reduction across various settings, including acoustic rooms, acoustic linings
<table><tr><td><b>5</b></td><td><b>Experimental validation</b></td><td><b>14</b></td></tr><tr><td>5.1</td><td>Experimental set-up</td><td>14</td></tr><tr><td>5.2</td><td>Relative displacement</td><td>15</td></tr><tr><td>5.3</td><td>Experimental and analytical results</td><td>16</td></tr></table>

<table><tr><td><b>6</b></td><td><b>Conclusions</b></td><td><b>17</b></td></tr></table>

<table><tr><td><b>References</b></td><td><b>18</b></td></tr></table>

## Nomenclature

Primary quantities introduced in the paper are recalled below with emphasis on the main explicit dependencies to other quantities. They are no repeated in the text. When deemed appropriate, secondary quantities with appropriate dependencies are directly introduced in the text.

### Geometric constants

<table><thead><tr><td>Symbol</td><td>Description</td><td>Explicit dependencies</td><td>Unit</td></tr></thead><tbody><tr><td>d</td><td>diameter of perforation</td><td>–</td><td>m</td></tr><tr><td>h</td><td>plate thickness</td><td>–</td><td>m</td></tr><tr><td>L<sub>x</sub></td><td>plate length in x direction</td><td>–</td><td>m</td></tr><tr><td>L<sub>y</sub></td><td>plate length in y direction</td><td>–</td><td>m</td></tr></tbody></table>

### Macroscopic parameters

<table><thead><tr><td>Symbol</td><td>Description</td><td>Explicit dependencies</td><td>Unit</td></tr></thead><tbody><tr><td>α<sub>∞</sub></td><td>effective tortuosity</td><td>φ, d</td><td>–</td></tr><tr><td>ε</td><td>end-length correction</td><td>φ, d, h</td><td>s m<sup>-1</sup></td></tr><tr><td>φ</td><td>perforation ratio</td><td>d</td><td>–</td></tr><tr><td>σ</td><td>nonlinear fluid-flow resistivity</td><td>φ, d, x, y, t</td><td>N s m<sup>-4</sup></td></tr><tr><td>σ<sub>0</sub></td><td>static fluid-flow resistivity</td><td>φ, d</td><td>N s m<sup>-4</sup></td></tr><tr><td>ξ</td><td>Forchheimer's parameter</td><td>φ, d, h</td><td>s m<sup>-1</sup></td></tr><tr><td>C<sub>d</sub></td><td>discharge coefficient</td><td>φ, d, h</td><td>–</td></tr></tbody></table>

### Mechanical parameters

<table><thead><tr><td>Symbol</td><td>Description</td><td>Explicit dependencies</td><td>Unit</td></tr></thead><tbody><tr><td>α</td><td>Biot's elastic coupling coefficient</td><td>–</td><td>–</td></tr><tr><td>α<sub>ac</sub></td><td>acoustic absorption</td><td>φ, d</td><td>–</td></tr><tr><td>η</td><td>loss factor</td><td>–</td><td>–</td></tr><tr><td>η<sub>i</sub></td><td>nonlinear loss factor of mode i</td><td>φ, d, ω, ξ, ω̇</td><td>–</td></tr><tr><td>η<sub>i</sub><sup>L</sup></td><td>linear loss factor of mode i</td><td>φ, d, ω</td><td>–</td></tr><tr><td>μ<sub>f</sub></td><td>fluid dynamic viscosity</td><td>–</td><td>kg m<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>ν</td><td>Poisson's coefficient</td><td>–</td><td>–</td></tr><tr><td>ρ</td><td>plate effective density</td><td>φ</td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>f</sub></td><td>fluid density</td><td>–</td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>s</sub></td><td>structure density</td><td>–</td><td>kg m<sup>-3</sup></td></tr><tr><td>D</td><td>bending stiffness</td><td>φ</td><td>GPa m<sup>3</sup></td></tr><tr><td>E</td><td>Young's modulus</td><td>–</td><td>GPa</td></tr><tr><td>K<sub>f</sub></td><td>fluid bulk modulus</td><td>–</td><td>kPa</td></tr><tr><td>M<sub>f</sub></td><td>Biot's elastic coefficient</td><td>–</td><td>N m<sup>-2</sup></td></tr><tr><td>R<sub>e</sub></td><td>Reynold's number</td><td>φ, d, ω̇</td><td>–</td></tr></tbody></table>

### Other parameters

<table><thead><tr><td>Symbol</td><td>Description</td><td>Explicit dependencies</td><td>Unit</td></tr></thead><tbody><tr><td>(x, y)</td><td>coordinate system attached to plate</td><td>–</td><td>m</td></tr><tr><td>(x<sub>e</sub>, y<sub>e</sub>)</td><td>coordinate of excitation point</td><td>–</td><td>m</td></tr><tr><td>(x<sub>g</sub>, y<sub>g</sub>)</td><td>coordinate of observation point</td><td>–</td><td>m</td></tr><tr><td>x</td><td>vector of Fourier coefficients</td><td>ω</td><td>m</td></tr><tr><td>ω</td><td>driving angular frequency</td><td>–</td><td>rad s<sup>-1</sup></td></tr><tr><td>ω<sub>c</sub></td><td>characteristic angular frequency</td><td>φ, d</td><td>rad s<sup>-1</sup></td></tr><tr><td>ω<sub>c</sub><sup>NL</sup></td><td>nonlinear characteristic angular frequency</td><td>φ, d</td><td>rad s<sup>-1</sup></td></tr><tr><td>ω<sub>i</sub></td><td>MPP angular frequency of mode i</td><td>–</td><td>rad s<sup>-1</sup></td></tr><tr><td>Ψ<sub>i</sub></td><td>modal shape of mode i</td><td>x, y</td><td>–</td></tr><tr><td>f<sub>NL</sub></td><td>nonlinear force in frequency domain</td><td>x, ω</td><td>N</td></tr><tr><td>f̃<sub>ext</sub></td><td>forcing force in frequency domain</td><td>ω</td><td>N</td></tr><tr><td>f<sub>NL</sub></td><td>nonlinear force in time domain</td><td>x, y, t</td><td>N</td></tr><tr><td>f<sub>c</sub></td><td>characteristic frequency</td><td>φ, d</td><td>Hz</td></tr><tr><td>F<sub>ext</sub></td><td>forcing force amplitude</td><td>–</td><td>N</td></tr><tr><td>f<sub>ext</sub></td><td>forcing force in time domain</td><td>x, y, t, ω</td><td>N</td></tr><tr><td>N</td><td>number of degrees-of-freedom</td><td>–</td><td>–</td></tr><tr><td>N<sub>h</sub></td><td>number of harmonics</td><td>–</td><td>–</td></tr></tbody></table>
within flow ducts [1], and even in nuclear engines and reactors [2]. These simple structures can also be
crafted from various materials, ensuring greater safety and environmental friendliness with reduced risk of
degradation. Moreover, these structures are well-suited for use in challenging environments characterized
by substantial excitations or high temperatures, such as in engine enclosures as acoustic barriers and
within sound suppressors, to provide a few examples.

Theoretical and experimental studies were carried out to characterize the linear response of the acoustic absorption of microperforated plates. Maa [3] proposed to model an MPP using Kirchhoff's equations. Atalla and Sgard [4] used the Johnson-Allard approach to model the microperforated plate as an equivalent fluid. They proposed a modified tortuosity to account for the distortion of fluid flow at the end of the perforations and the interaction between them. The observed acoustic absorption was induced by dissipative frictional mechanisms in the thermal and viscous boundary layers for a perforation diameter of the order of their thickness. In the area of structural dynamics, Gallerand et al. [5, 6] proposed an analytical model, validated by experimental measurements, to explore the dynamic response of a finite-size MPP saturated by air. It was shown that additional damping, occurs around a characteristic frequency, because of dissipative mechanisms appearing in the thermoviscous boundary layers at the fluid-solid interface of the perforations.

However, microperforated plates are likely to be used in harsh environments subject to strong mechanical and/or acoustic excitations. This is why MPP, traditionally used as acoustic materials, have been studied under high levels of acoustic excitation. In the context of nonlinear acoustics, the acoustic velocity of particles in perforations can attain high values. This results in flow separation and vortices at the sharp edges of the perforations. The formation of vortices absorbs the acoustic wave energy and thus increases the acoustic resistance of the MPP's perforation [7, 8, 9]. Acoustic impedance is then modified to accommodate nonlinear effects. In this context, Maa [9] proposed a formulation of the acoustic impedance for an MPP at a high pressure level and found a relationship between the nonlinear resistance and the Mach number defined from the particle velocity in the perforation. Cummings [10, 11] proposed a model to predict the temporal behavior of reactance (the imaginary part of the acoustic impedance). It was shown that the energy absorbed by the material depends on the sound level but that there is no linear relationship between these two parameters. Melling [12] provided a nonlinear reactance model based on a quasi-stationary approximation of acoustic flow through perforations. Experimentally, Ingard and Ising [13] showed that, in the nonlinear regime, the fluid separates at the perforation outlet in the form of jets. They also showed that pressure and velocity amplitudes were linked by a square-law relationship. Forchheimer [14] proposed an empirical relationship between velocity and pressure to correct for resistive effects when fluid flow rates in porous media become significant. The relationship is an extension of Darcy's linear law for large fluid flow rates in porous media. At the viscous term of Darcy's law, this expression includes a quadratic or cubic inertial term [14, 15]. Adler et al. [16] showed that there is no quadratic correction to Darcy's law for small Reynolds numbers $R_e$ (less than unity). Conversely, quadratic correction dominates in the case where $R_e > 1$. In the context of MPP, $R_e > 1$ is expected for airflow velocities greater than $1 \text{ mm s}^{-1}$ in millimeter perforations.

The law proposed by Forchheimer has been applied and used for MPP under high-acoustic excitations. Using the equivalent fluid model and the approach proposed by [4], this previously defined approach for porous materials can be applied to MPP. The Forchheimer correction was used by Tayong et al. [7] to investigate the acoustic absorption of MPP supported by an air cavity. They demonstrated that absorption coefficient increases with the Mach number up to a critical Mach number. If the critical Mach number of an MPP is too low (meaning that the critical Mach number is within the linear amplitude regime), the absorption coefficient of the MPP only decreases with the Mach number in the nonlinear amplitude regime. In another work, Tayong et al. [17] experimentally explored the interaction between multiple perforations in a nonlinear acoustic regime. Forchheimer's quadratic law was also used by Laly et al. [18] who proposed a model to characterize the acoustic impedance of an MPP at high acoustic pressure. They modeled the microperforated plate as an equivalent fluid using the Johnson-Allard approach.

In all of the works presented above, the MPP is considered rigid. For low excitation amplitudes (acoustic and mechanical), there are approaches that consider structural vibrations of the MPP. However, no studies have examined the combined vibrational and acoustic behavior of a flexible MPP in the nonlinear excitation regime. In order to grasp the nonlinearity domain of the present study, it is nevertheless important to emphasize that numerous studies are available on the response of conventional plates or perforated plates [19]. These studies assume that there is no interaction between the fluid in the perforations and the structure, subjected to high amplitudes. In linear theories, it is assumed that the deformations are small compared to the thickness and that the shear deformations in the straight section can be neglected. This assumption is no longer valid for large-amplitude excitations. Consequently, the derivation of governing equations with large deflections requires special attention. Von Kármán's nonlinear plate theory was widely used to explore the forced harmonic vibration of plates, beams at large amplitudes [20, 21, 22, 23] or more complicated geometry like a plate with acoustic black hole [24, 25].
<table><tr><td>t</td><td>time</td><td>-</td><td>s</td></tr><tr><td>V</td><td>amplitude of ṇ</td><td>x, y, φ, d</td><td>m s<sup>-1</sup></td></tr><tr><td>V<sub>cr</sub></td><td>critical relative velocity</td><td>φ, d, ω</td><td>m s<sup>-1</sup></td></tr><tr><td>w</td><td>fluid-solid relative displacement</td><td>x, y, t</td><td>m</td></tr><tr><td>w<sub>s</sub></td><td>solid displacement</td><td>x, y, t</td><td>m</td></tr></table>

# 1. Introduction

Microperforated plates (MPP) are commonly used as lightweight acoustic-absorbing materials, serving as a viable substitute for porous materials and conventional acoustic resonators. MPP have versatile applications in acoustic reduction across various settings, including acoustic rooms, acoustic linings within flow ducts [1], and even in nuclear engines and reactors [2]. These simple structures can also be crafted from various materials, ensuring greater safety and environmental friendliness with reduced risk of degradation. Moreover, these structures are well suited for use in challenging environments characterized by substantial excitations or high temperatures, such as in engine enclosures as acoustic barriers and within sound suppressors, to provide a few examples.

Theoretical and experimental studies were carried out to characterize the linear response of the acoustic absorption of microperforated plates. Maa [3] proposed to model an MPP using Kirchhoff's equations. Atalla and Sgard [4] used the Johnson-Allard approach to model the microperforated plate as an equivalent fluid. They proposed a modified tortuosity to account for the distortion of fluid flow at the end of the perforations and the interaction between them. The observed acoustic absorption was induced by dissipative frictional mechanisms in the thermal and viscous boundary layers for a perforation diameter of the order of their thickness. In the area of structural dynamics, Gallerand et al. [5, 6] proposed an analytical model, validated by experimental measurements, to explore the dynamic response of a finite-size MPP saturated by air. It was shown that additional damping, occurs around a characteristic frequency, because of dissipative mechanisms appearing in the thermoviscous boundary layers at the fluid-solid interface of the perforations. In summary, MPP can be used in two different ways: (1) coupling MPP to air cavities, where the acoustic resonance between the MPP and the cavity provides the sound absorption, and (2) microperforating an existing structure to add structural damping on modes.

However, microperforated plates are likely to be used in harsh environments subject to strong mechanical and/or acoustic excitations. This is why MPP, traditionally used as acoustic materials, have been studied under high levels of acoustic excitation. In the context of nonlinear acoustics, the acoustic velocity of particles in perforations can attain high values. This results in flow separation and vortices at the sharp edges of the perforations. The formation of vortices absorbs the acoustic wave energy and thus increases the acoustic resistance of the MPP's perforation [7-9]. The acoustic impedance is then modified to accommodate nonlinear effects. In this context, Maa [9] proposed a formulation of the acoustic impedance for an MPP at a high-pressure level and found a relationship between the nonlinear resistance and the Mach number defined from the particle velocity in the perforation. Cummings [10, 11] proposed a model to predict the temporal behavior of reactance (the imaginary part of the acoustic impedance). It was shown that the energy absorbed by the material depends on the sound level but that there is no linear relationship between these two parameters. Melling [12] provided a nonlinear reactance model based on a quasi-stationary approximation of acoustic flow through perforations. Experimentally, Ingard and Ising [13] showed that, in the nonlinear regime, the fluid separates at the perforation outlet in the form of jets. They also showed that pressure and velocity amplitudes were linked by a square-law relationship. Forchheimer [14] proposed an empirical relationship between velocity and pressure to correct for resistive effects when fluid flow rates in porous media become significant. The relationship is an extension of Darcy's linear law for large fluid flow rates in porous media. At the viscous term of Darcy's law, this expression includes a quadratic or cubic inertial term [14, 15]. Adler et al. [16] showed that there is no quadratic correction to Darcy's law for small Reynolds numbers<fnref n="1" /> $R_e$ (less than unity). Conversely, quadratic correction dominates in the case where $R_e > 1$. In the context of MPP, $R_e > 1$ is expected for airflow velocities greater than 1 mm s⁻<fnref n="1" /> in millimeter perforations.

The law proposed by Forchheimer has been applied and used for MPP under high-acoustic excitations. Using the equivalent fluid model and the approach proposed by [4], this previously defined approach for porous materials can be applied to MPP. The Forchheimer correction was used by Tayong et al. [7] to investigate the acoustic absorption of MPP supported by an air cavity. They demonstrated that the absorption coefficient increases with the Mach number up to a critical Mach number. If the critical Mach number of an MPP is too low (meaning that the critical Mach number is within the linear amplitude

<footnote n="1">¹The Reynolds number as well as the Mach number are based on the acoustic wave. The fluid flow velocity used to determine these addimentioned parameters is the root-mean-square velocity, not the instantaneous fluid flow velocity.</footnote>
In these studies, only structural/geometric nonlinearities are considered. For MPP saturated with air in a dynamic context, we therefore distinguish two types of nonlinearity: (i) acoustic nonlinearity, and (ii) geometric nonlinearity. The geometric nonlinearity can be captured by cubic stiffness terms [22, 23], while acoustic nonlinearity can be captured by quadratic stiffness terms [7]. As the magnitude of the fluid particle velocity increases in the perforation, the acoustic nonlinearity may be activated first, before the geometric nonlinearity. Indeed, via viscous fluid-structure coupling, the solid will impose its motion on the fluid. This coupling is strongly affected by the thickness of the visco-thermal boundary layers. Thus, for perforation diameters largely smaller than the thickness of the visco-thermal boundary layers, the fluid is totally “pulled” out of the perforation. Its velocity in the perforation is therefore zero. For perforation diameters significantly greater than the thickness of the boundary layers, fluid movement in the center of the perforation will not be affected by solid movement. In some cases, the fluid velocity is also zero. In these cases, the relative velocity between the fluid and the solid is equal to the velocity of the solid, and the acoustic nonlinearity linked to the fluid-structure coupling will not appear. Between these two limits, the reduction in air flow surface during the movement of the MPP results in a fluid velocity in the perforations greater than that of the solid. Acoustic particle velocity increases significantly when the wave is compressed in the microperforations [13]. Nonlinearities caused by fluid-structure coupling (related to the fluid) will thus appear before those related to the solid. In this paper, we are only interested in the acoustic nonlinearity, also called Forchheimer’s nonlinearity, due to the fluid into the perforations. In the Forchheimer regime, a high harmonic flow velocity through the perforations is assumed.

This article is organized into four sections. First, Section 2 proposes a nonlinear vibratory model based on the linear vibratory MPP model adapted to a nonlinear acoustic framework using Forchheimer’s law. This leads to a system of two coupled partial differential equations (PDE). Secondly, Section 3 provides details on the procedure for solving the resulting nonlinear PDE. Thirdly, Section 4 analyzes the effects of the nonlinear damping on the MPP behaviour. Finally, Section 5 presents experimental results to validate the developments.

## 2. Model

This section recalls the equations governing the dynamics of a finite size microperforated plate saturated by a light fluid, as proposed in [5]. They are then extended to the nonlinear framework where the velocity of the fluid in the microperforations is large and the inertial effects in the microperforations are no longer negligible. The latter are captured by a correction term proposed by Forchheimer [14, 7]. This correction leads to an additional antisymmetric quadratic damping term in the governing equations.

### 2.1. Dynamics of a finite-size microperforated plate saturated by a light fluid

A simply supported finite-size MPP of dimension $L_x \times L_y \times h$ placed in the $xy$ plane is investigated. Only displacements along the $z$-axis is considered. The plate is excited by an external force $f_{ext}$. The equations of motion read [5]

$$h(\rho\ddot{w}_s + \rho_f\ddot{w}) + \left(D + \frac{\alpha^2 M_f h^3}{12}\right)\nabla^4 w_s = f_{ext} \quad (1a)$$

$$\rho_f \ddot{w}_s + \frac{\rho_f \alpha_\infty}{\phi} \ddot{w} + \sigma_0 \dot{w} + \alpha M_f \nabla^2 w_s = 0, \quad (1b)$$

where $\nabla$ is the classical gradient differential operator. They are obtained by considering the MPP as a particular porous plate [4] using an alternative form of Biot’s theory [26]. Equation (1) governs the vibratory response of the MPP in the linear regime where $w_s$ is the solid motion and $w = \phi(w_s - w_f)$ is the relative fluid-solid motion with $\phi$ the perforation ratio and $w_f$ the fluid motion. Equation (1a) represents the elastic response of the homogeneous solid plate, while Equation (1b) corresponds to the relative fluid-solid motion. The fluid-solid density mixture $\rho = (1 - \phi)\rho_s + \phi\rho_f$ is a function of the solid density $\rho_s$ and fluid density $\rho_f$. The parameters $\alpha$ and $M_f$ are elastic coefficients defined by Biot [27]. In the context of an MPP saturated by a lightweight fluid the Biot’s coefficients are simply written $\alpha = \phi$ and $M_f = K_f/\phi$ [26] where $K_f$ is the bulk modulus of air. The coefficient $D$ is the bending stiffness. In order to consider the influence of the microperforations in the MPP stiffness, $D$ is updated to [28, 29]

$$D = \frac{ECh^3}{12(1 - \nu^2)} \quad \text{with} \quad C = \frac{(1 - \phi)^2}{1 + (2 - 3\nu)\phi}, \quad (2)$$

where $E$ and $\nu$ are respectively the Young’s modulus and Poisson’s coefficient for the non-perforated plate. For a vibrating MPP saturated with a lightweight fluid, the elastic modulus of the plate is much larger than that of the fluid, i.e. $E \gg M_f$. In this case, with MPP vibrating in air, this condition is fulfilled and only inertial interactions are accounted for in the plate equations. In Equation (1), the stiffness coefficient,
regime), the absorption coefficient of the MPP only decreases with the Mach number in the nonlinear amplitude regime. In another work, Tayong et al. [17] experimentally explored the interaction between multiple perforations in a nonlinear acoustic regime. Forchheimer's quadratic law was also used by Laly et al. [18] who proposed a model to characterize the acoustic impedance of an MPP at high acoustic pressure. They modeled the microperforated plate as an equivalent fluid using the Johnson-Allard approach.

In all of the works presented above, the MPP is considered rigid. For low excitation amplitudes (acoustic and mechanical), there are approaches that consider structural vibrations of the MPP. However, no studies have examined the combined vibrational and acoustic behavior of a flexible MPP in the nonlinear excitation regime. In order to grasp the nonlinearity domain of the present study, it is nevertheless important to emphasize that numerous studies are available on the response of conventional plates or perforated plates [19]. These studies assume that there is no interaction between the fluid in the perforations and the structure, subjected to high amplitudes. In linear theories, it is assumed that the deformations are small compared to the thickness and that the shear deformations in the straight section can be neglected. This assumption is no longer valid for large-amplitude excitations. Consequently, the derivation of governing equations with large deflections requires special attention. Von Kármán's nonlinear plate theory was widely used to explore the forced harmonic vibration of plates, beams at large amplitudes [20-23] or more complicated geometry like a plate with acoustic black hole [24, 25].

In these studies, only structural/geometric nonlinearities are considered. For MPP saturated with air in a dynamic context, we therefore distinguish two types of nonlinearity: (i) acoustic nonlinearity, and (ii) geometric nonlinearity. The geometric nonlinearity can be captured by cubic stiffness terms [22, 23], while acoustic nonlinearity can be captured by quadratic velocity terms [7]. As the magnitude of the fluid particle velocity increases in the perforation, the acoustic nonlinearity may be activated first, before the geometric nonlinearity. Indeed, via viscous fluid-structure coupling, the solid will impose motion on the fluid. This coupling is strongly affected by the thickness of the visco-thermal boundary layers. Thus, for perforation diameters largely smaller than the thickness of the visco-thermal boundary layers, the fluid is totally “pulled” out of the perforation. Its velocity in the perforation is therefore zero. For perforation diameters significantly greater than the thickness of the boundary layers, fluid movement in the center of the perforation will not be affected by solid movement. In some cases, the fluid velocity is also zero. In these cases, the relative velocity between the fluid and the solid is equal to the velocity of the solid, and the acoustic nonlinearity linked to the fluid-structure coupling will not appear. Between these two limits, the reduction in air flow surface during the movement of the MPP results in a fluid velocity in the perforations greater than that of the solid. Acoustic particle velocity increases significantly when the wave is compressed in the microperforations [13]. Nonlinearities caused by fluid-structure coupling (related to the fluid) will thus appear before those related to the solid. In this paper, we are only interested in the acoustic nonlinearity, also called Forchheimer's nonlinearity, due to the fluid into the perforations. In the Forchheimer regime, a high harmonic flow velocity through the perforations is assumed.

This article is organized into four sections. First, Section 2 proposes a nonlinear vibratory model based on the linear vibratory MPP model adapted to a nonlinear acoustic framework using Forchheimer's law. This leads to a system of two coupled partial differential equations (PDE). Secondly, Section 3 provides details on the procedure for solving the resulting nonlinear PDE. Thirdly, Section 4 analyzes the effects of the nonlinear damping on the MPP behaviour. Finally, Section 5 presents experimental results to validate the developments.

## 2. Model

This section recalls the equations governing the dynamics of a finite size microperforated plate saturated by a light fluid, as proposed in [5]. They are then extended to the nonlinear framework where the velocity of the fluid in the microperforations is large and the inertial effects in the microperforations are no longer negligible. The latter are captured by a correction term proposed by Forchheimer [7, 14]. This correction leads to an additional antisymmetric quadratic damping term in the governing equations.

### 2.1. Dynamics of a finite-size microperforated plate saturated by a light fluid

A simply supported finite-size MPP of dimensions $L_x \times L_y \times h$ placed in the $xy$ plane is investigated. Only displacements along the $z$-axis is considered. The plate is excited by an external force $f_{ext}$. The equations of motion read [5]

$$h(\rho\ddot{w}_s + \rho_f\ddot{w}) + \left(D + \frac{\alpha^2 M_f h^3}{12}\right)\nabla^4 w_s = f_{ext} \quad (1a)$$

$$\rho_f \ddot{w}_s + \frac{\rho_f \alpha \infty}{\phi} \ddot{w} + \sigma_0 \dot{w} + \alpha M_f \nabla^2 w_s = 0 \quad (1b)$$
$D + \frac{\alpha^2 M_f h^3}{12}$ can be approximated by $D$, as considered in the rest of the present paper. Also, the bulk modulus of air $K_f$ is considered constant and equal to its adiabatic value in the resolution of the system of Equation (1). All Johnson-Champoux-Allard (JCA) parameters defined for a porous medium can be rewritten for an MPP as a function of $\phi$ and the diameter of the perforations $d$. Thus, the static airflow resistivity and the tortuosity are defined by

$$ \sigma_0 = \frac{32\mu_f}{\phi d^2} \quad \text{and} \quad \alpha_\infty = 1 + \frac{2\epsilon}{h}, \qquad (3) $$

where $\mu_f$ is the fluid dynamic viscosity and $\epsilon = 0.24\sqrt{\pi d^2}(1 - 1.14\sqrt{\phi})$ [4] is an end correction factor used to consider the fluid radiation inside the perforations and the distortion of the fluid flow at the perforation orifices.

## 2.2. Added damping in linear regime

In a linear context, an MPP exhibits viscous added damping due to energy dissipation induced by fluid-structure interactions coupled with thermo-viscous ones between shearing adjacent fluid layers near the perforation solid walls. The loss factor resulting of these mechanisms writes in terms of $\omega$, the angular frequency [5],

$$ \eta_i^L = - \frac{h\omega I_i^{(2)} I_i^{(3)}}{D I_i^{(1)}} \frac{\alpha M_f \phi^2 \rho_f \sigma_0 I_i^{(3)}}{(\alpha_\infty \rho_f \omega I_i^{(3)})^2 + (\phi \sigma_0 I_i^{(3)})^2}, \qquad (4) $$

where $I_i^{(p)}$ with $p = 1, 2, 3$ is the $i$th component of the projection matrix $\mathbf{I}^{(k)}$ defined from $\Psi(x, y)$ which stores the $\Psi_i(x, y)$ term corresponding to the $i$th mode shape function such that<fnref n="1" />

$$ \mathbf{I}^{(1)} = \iint_S \nabla^2 \Psi(x, y) \nabla^2 \Psi^T(x, y) dxdy, \quad \mathbf{I}^{(2)} = \iint_S \nabla^2 \Psi(x, y) \Psi^T(x, y) dxdy, \quad (5) \\ \mathbf{I}^{(3)} = \iint_S \Psi(x, y) \Psi^T(x, y) dxdy. $$

The loss factor, provided in Equation (4), reaches a maximum at a characteristic angular frequency $\omega_c$ given by [5]

$$ \omega_c = \frac{32\mu_f}{\alpha_\infty \rho_f d^2} \quad \text{and its corresponding frequency} \quad f_c = \frac{16\mu_f}{\pi\alpha_\infty \rho_f d^2}, \qquad (6) $$

defined from Biot's frequency for porous materials [30, 26] and adapted for the MPP case. The added damping delivered by the microperforations is maximum at the characteristic frequency, but also has an effect in a range of frequencies around $f_c$. Since $f_c$ is a function of $d$, the perforation diameter can be adapted to induce maximum added damping at a resonance frequency of the MPP, i.e. to make $f_c$ coincide with a natural frequency of the MPP.

## 2.3. Added damping in nonlinear regime

When the fluid velocity in the microperforations becomes sufficiently high, the airflow resistive effects occurring in the microperforations are no longer constant and increase linearly with the fluid velocity [7, 31]. They are incorporated into the model using Forchheimer's law, which is defined for a rigid porous material subjected to high flow velocity in the pores. According to Forchheimer's law, the airflow resistivity reads [15, 18]

$$ \sigma = (1 + \varepsilon|\dot{w}_f|)\sigma_0. \qquad (7) $$

with $\varepsilon$, the Forchheimer's parameter. In the context of a vibrating MPP, the relative fluid-solid velocity corresponds to the fluid velocity for a rigid frame MPP, that is $w$ corresponds to $w_f$, and Equation (7) becomes

$$ \sigma = (1 + \varepsilon|\dot{w}|)\sigma_0. \qquad (8) $$

In fact, the absolute fluid velocity used in Equation (7) corresponds to the relative velocity between fluid and solid in the case of a vibrating MPP.

<footnote n="1">¹More details about the spatial and temporal discretization of the PDE are provided in Section 3.</footnote>
where $\nabla$ is the classical gradient differential operator. They are obtained by considering the MPP as a particular porous plate [4] using an alternative form of Biot's theory [26]. Equation (1) governs the vibratory response of the MPP in the linear regime where $w_s$ is the solid motion and $w = \phi(w_s - w_f)$ is the relative fluid-solid motion with $\phi$ the perforation ratio and $w_f$ the fluid motion. Equation (1a) represents the elastic response of the homogeneous solid plate, while Equation (1b) corresponds to the relative fluid-solid motion. The fluid-solid density mixture $\rho = (1 - \phi)\rho_s + \phi\rho_f$ is a function of the solid density $\rho_s$ and fluid density $\rho_f$. The parameters $\alpha$ and $M_f$ are elastic coefficients defined by Biot [27]. In the context of an MPP saturated by a lightweight fluid the Biot's coefficients are simply written $\alpha = \phi$ and $M_f = K_f/\phi$ [26] where $K_f$ is the bulk modulus of air. The coefficient $D$ is the bending stiffness. In order to consider the influence of the microperforations in the MPP stiffness, $D$ is updated to [28, 29]

$$D = \frac{ECh^3}{12(1 - \nu^2)} \quad \text{with} \quad C = \frac{(1 - \phi)^2}{1 + (2 - 3\nu)\phi}, \qquad (2)$$

where $E$ and $\nu$ are respectively the Young's modulus and Poisson's coefficient for the non-perforated plate. For a vibrating MPP saturated with a lightweight fluid, the elastic modulus of the plate is much larger than that of the fluid, i.e. $E \gg M_f$. In this case, with MPP vibrating in air, this condition is fulfilled and only inertial interactions are accounted for in the plate equations. In Equation (1), the stiffness coefficient, $D + \frac{\alpha^2 M_f h^3}{12}$ can be approximated by $D$, as considered in the rest of the present paper. Also, the bulk modulus of air $K_f$ is considered constant and equal to its adiabatic value in the resolution of the system of Equation (1). All Johnson-Champoux-Allard (JCA) parameters defined for a porous medium can be rewritten for an MPP as a function of $\phi$ and the diameter of the perforations $d$. In this paper, all perforations are assumed to be cylindrical with the same circular cross-section, the size of which can be varied depending on the investigated case. The static airflow resistivity and the tortuosity are defined by

$$\sigma_0 = \frac{32\mu_f}{\phi d^2} \quad \text{and} \quad \alpha_\infty = 1 + \frac{2\epsilon}{h} \qquad (3)$$

where $\mu_f$ is the fluid dynamic viscosity and $\epsilon = 0.24\sqrt{\pi d^2}(1 - 1.14\sqrt{\phi})$ [4] is an end correction factor used to incorporate the fluid radiation inside the perforations and the distortion of the fluid flow at the perforation orifices.

## 2.2. Added damping in linear regime

In a linear context, an MPP exhibits viscous added damping due to energy dissipation induced by fluid-structure interactions coupled with thermo-viscous ones between shearing adjacent fluid layers near the perforation solid walls. The loss factor resulting of these mechanisms is written in terms of the angular frequency $\omega$ [5],

$$\eta_i^L(\omega) = -\frac{h\omega I_i^{(2)} I_i^{(3)}}{D I_i^{(1)}} \frac{\alpha M_f \phi^2 \rho_f \sigma_0 I_i^{(3)}}{(\alpha_\infty \rho_f \omega I_i^{(3)})^2 + (\phi \sigma_0 I_i^{(3)})^2}, \qquad (4)$$

where $I_i^{(p)}$ with $p = 1, 2, 3$ is the $i$th component of the projection matrix $\mathbf{I}^{(k)}$ defined from $\Psi(x, y)$ which stores the $\Psi_i(x, y)$ term corresponding to the $i$th mode shape function such that$^2$

$$\mathbf{I}^{(1)} = \iint_S \nabla^2 \Psi(x, y) \nabla^2 \Psi^T(x, y) \, dx \, dy, \quad \mathbf{I}^{(2)} = \iint_S \nabla^2 \Psi(x, y) \Psi^T(x, y) \, dx \, dy, \\ \mathbf{I}^{(3)} = \iint_S \Psi(x, y) \Psi^T(x, y) \, dx \, dy. \qquad (5)$$

The loss factor in Equation (4) reaches a maximum at the characteristic angular frequency [5]

$$\omega_c = \frac{32\mu_f}{\alpha_\infty \rho_f d^2} \quad \text{and its corresponding frequency } f_c = \frac{16\mu_f}{\pi\alpha_\infty \rho_f d^2}, \qquad (6)$$

defined from Biot's frequency for porous materials [26, 30] and adapted to the MPP case. The added damping delivered by the microperforations is maximum at the characteristic frequency, but also has an effect in a range of frequencies around $f_c$. Since $f_c$ is a function of $d$, the perforation diameter can be adapted to induce maximum added damping at a resonance frequency of the MPP, i.e. to make $f_c$ coincide with a natural frequency of the MPP.

<footnote>$^2$More details about the spatial and temporal discretization of the PDE are provided in Section 3 and in [5].</footnote>
The Forchheimer parameter is obtained through airflow resistivity measurement on a high-flow resistivity meter (Mecanum Inc. HF-SIGMA high-flow resistivity-meter). The Forchheimer parameter is determined by measuring the pressure drop along the sample at several flow rates and using straight-line interpolation at zero flow velocity, in accordance with Equation (7) [32, 33]. Through this procedure, measurements on an MPP were conducted with $d = 2.8$ mm, $\phi = 10\%$ and $h = 1$ mm yield$^2 \varepsilon = 1.43 \pm 0.025$ s m$^{-1}$. Another way to determine $\varepsilon$ is to use the empirical expression [18]

$$ \varepsilon = \frac{\beta \rho_f (1 - \phi^2)}{\pi h \phi C_d^2 \sigma_0} \phi, \qquad (9) $$

where $C_d$ is the discharge coefficient, commonly with $C_d \in 0.6 - 0.8$ and set here to 0.76. The unitless constant $\beta$ is assumed to be equal to 1.6 according to [18]. Using Equation (9) to determine $\varepsilon$ for an MPP $d = 2.8$ mm, $\phi = 10\%$ and $h = 1$ mm provide $\varepsilon = 1.42$ s m$^{-1}$. Comparison between Equation (9) and the experimental method suggests that Equation (9) appears a reasonable estimate for $\varepsilon$. Note that the Forchheimer's parameter depends only on the perforation parameters ($d$ and $\phi$) and the MPP thickness.

From Forchheimer's law in Equation (1) is derived the system

$$ h(\rho\ddot{w}_s + \rho_f\ddot{w}) + D\nabla^4 w_s = f_{\text{ext}}, \qquad (10a) $$

$$ \rho_f \ddot{w}_s + \frac{\rho_f \alpha_\infty}{\phi} \ddot{w} + \sigma_0 \dot{w} + \sigma_0 \varepsilon |\dot{w}| \dot{w} + \alpha M_f \nabla^2 w_s = 0, \qquad (10b) $$

in which temporal and spatial dependencies were dropped for readability. A supplementary quadratic damping term $\sigma_0 \varepsilon |\dot{w}| \dot{w}$, equivalent to a modal force always opposite to the direction of the velocity and proportional to its square, is incorporated.

### 3. Resolution procedure

#### 3.1. Harmonic method balance

This subsection provides details on the discretization procedure and the resolution of Equation (10). Solutions for $w_s$, $w$ and $f_{\text{ext}}$ are developed on the basis of orthogonal or quasi-orthogonal eigenfunctions $\Psi$. For instance, the plate displacement is assumed to be of the form

$$ w_s(x, y, t) = \sum_{i=1}^{N} w_i^s(t) \Psi_i(x, y), \qquad (11) $$

where $w_i^s$ is the generalized coordinate of modeshape $\Psi_i$, $i \in \mathbb{N}^*$ and $N$ is the number of degrees-of-freedom (dof). Index $i$ corresponds to a reorganization of the pair of eigenmodes $(m, n) \in (\mathbb{N}^*, \mathbb{N}^*)$. Similar expressions are assumed for $w$ and $f_{\text{ext}}$ is chosen as a force applied at the point $(x_e, y_e)$ and is written for the $i$th mode such as

$$ f_{\text{ext}}(x_e, y_e, t) = F_{\text{ext}} \sum_{i=1}^{N} \Psi_i(x_e, y_e) q_i^{\text{ex}}(t). \qquad (12) $$

After projection on eigenmode basis of a non-perforated plate, the matrix form of the dynamics reads

$$ \mathbf{M}_1 \ddot{\mathbf{w}}_s(t) + \mathbf{M} \ddot{\mathbf{w}}(t) + \mathbf{K}_1 \mathbf{w}_s(t) = \mathbf{f}_{\text{ex}}(t), \qquad (13a) $$

$$ \mathbf{M} \ddot{\mathbf{w}}_s(t) + \mathbf{M}_2 \ddot{\mathbf{w}}(t) + \mathbf{C} \dot{\mathbf{w}}(t) + \mathbf{f}_{\text{NL}}(\dot{\mathbf{w}}(t)) + \mathbf{K}_2 \mathbf{w}_s(t) = \mathbf{0}, \qquad (13b) $$

where

$$ \mathbf{K}_1 = D \int_S \nabla^2 \mathbf{\Psi}(x, y) \nabla^2 \mathbf{\Psi}^\top(x, y) dx dy, \qquad \mathbf{K}_2 = h \alpha M_f \int_S \nabla^2 \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) dx dy, \\ \mathbf{C} = h \sigma_0 \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) dx dy, \qquad \mathbf{M} = h \rho_f \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) dx dy, \qquad (14) \\ \mathbf{M}_1 = h \rho \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) dx dy, \qquad \mathbf{M}_2 = \frac{h \rho_f \alpha_\infty}{\phi} \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) dx dy. $$

<footnote>$^2$The perforation parameters and plate thickness used to experimentally determine the Forchheimer parameter correspond to those used in the experimental validation section (Section 5).</footnote>
### 2.3. Added damping in nonlinear regime

When the fluid velocity in the microperforations becomes sufficiently high, the airflow resistive effects occurring in the microperforations are no longer constant and increase linearly with the fluid velocity [7, 31]. They are incorporated into the model using Forchheimer's law, which is defined for a rigid porous material subjected to high flow velocity in the pores. According to Forchheimer's law, the airflow resistivity reads [15, 18]

$$ \sigma = (1 + \xi|\dot{w}_f|)\sigma_0. \qquad (7) $$

with Forchheimer's parameter $\xi$. In the context of a vibrating MPP, the relative fluid-solid velocity corresponds to the fluid velocity for a rigid frame MPP, studied in an acoustic context, that is $w$ corresponds to $w_f$, and Equation (7) becomes

$$ \sigma = (1 + \xi|\dot{w}|)\sigma_0. \qquad (8) $$

In fact, the absolute fluid velocity in Equation (7) corresponds to the relative velocity between the fluid and solid in the case of a vibrating MPP.

The Forchheimer parameter is obtained through airflow resistivity measurement on a high-flow resistivity meter (Mecanum Inc. HF-SIGMA high-flow resistivity-meter). The Forchheimer parameter is determined by measuring the pressure drop along the sample at several flow rates and using straight-line interpolation at zero flow velocity, in accordance with Equation (7) [32, 33]. Through this procedure, measurements on an MPP were conducted with $d = 2.8$ mm, $\phi = 10\%$ and $h = 1$ mm yielding$^3 \xi = 1.43 \pm 0.025$ s m$^{-1}$. Another way to determine $\xi$ is to use the empirical expression [18]

$$ \xi = \frac{\beta \rho_f (1 - \phi^2)}{\pi h \phi C_d^2 \sigma_0} \phi, \qquad (9) $$

where $C_d$ is the discharge coefficient, commonly with $C_d \in 0.6 - 0.8$. A typical value is $C_d \approx 0.76$ from [18]. The unitless constant $\beta$ is assumed to be equal to 1.6 according to [18]. Equation (9) with the above values for $d$, $\phi$ and $h$ leads to $\xi = 1.42$ s m$^{-1}$, in good agreement with the experimental value. Note that the Forchheimer parameter depends only on the perforation parameters ($d$ and $\phi$) and the MPP thickness.

From Forchheimer's law in Equation (1) is derived the system

$$ h(\rho\ddot{w}_s + \rho_f\ddot{w}) + D\nabla^4 w_s = f_{\text{ext}}, \qquad (10a) $$

$$ \rho_f \ddot{w}_s + \frac{\rho_f \alpha_\infty}{\phi} \ddot{w} + \sigma_0 \dot{w} + \sigma_0 \xi |\dot{w}| \dot{w} + \alpha M_f \nabla^2 w_s = 0, \qquad (10b) $$

in which temporal and spatial dependencies were dropped for readability. A supplementary quadratic damping term $\sigma_0 \xi |\dot{w}| \dot{w}$, equivalent to a modal force always opposite to the direction of the velocity and proportional to its square, is incorporated.

## 3. Solution method

### 3.1. Harmonic balance method

This subsection provides details on the discretization procedure and the resolution of Equation (10). Solutions for $w_s$, $w$ and $f_{\text{ext}}$ are developed on the basis of orthogonal or quasi-orthogonal eigenfunctions $\Psi$. For instance, the plate displacement is assumed to be of the form

$$ w_s(x, y, t) = \sum_{i=1}^{N} w_{s_i}(t) \Psi_i(x, y), \qquad (11) $$

where $w_{s_i}$ is the generalized coordinate of mode shape $\Psi_i$, $i \in \mathbb{N}^*$ and $N$ is the number of degrees-of-freedom (dof). Index $i$ corresponds to a reorganization of the pair of eigenmodes $(m, n) \in (\mathbb{N}^*, \mathbb{N}^*)$. Similar expressions are assumed for $w$, and $f_{\text{ext}}$ is chosen as a force applied at the point $(x_e, y_e)$, expressed as

$$ f_{\text{ext}}(x_e, y_e, t) = F_{\text{ext}} \sum_{i=1}^{N} \Psi_i(x_e, y_e) q_i^{\text{ex}}(t). \qquad (12) $$

<footnote>$^3$The perforation parameters and plate thickness used to experimentally determine the Forchheimer parameter correspond to those used in the experimental validation section (Section 5).</footnote>
where $S = L_x \times L_y$. Entry $i$ of the nonlinear force vector $\mathbf{f}_{\text{NL}}(\dot{\mathbf{w}}(t))$ reads

$$f_i^{\text{NL}}(\dot{\mathbf{w}}(t)) = h\sigma_0\epsilon \int_S \sum_{j=1}^{N} \dot{w}_j(t)\Psi_j(x, y) \left| \sum_{p=1}^{N} \dot{w}_p(t)\Psi_p(x, y) \right| \Psi_i(x, y)dxdy. \quad (15)$$

The Harmonic Balance Method (HBM) is implemented and the contributions are assumed to be periodic and therefore expressed as Fourier series truncated at the harmonic $N_h$ [34, 35, 36]:

$$\mathbf{w}_s(t) = \frac{1}{2}\mathbf{A}_0 + \sum_{k=1}^{N_h} \mathbf{A}_k \cos k\omega t + \mathbf{B}_k \sin k\omega t \quad \text{and} \quad \dot{\mathbf{w}}(t) = \frac{1}{2}\mathbf{C}_0 + \sum_{k=1}^{N_h} \mathbf{C}_k \cos k\omega t + \mathbf{D}_k \sin k\omega t. \quad (16)$$

The vectors $\mathbf{A}_k, \mathbf{B}_k, \mathbf{C}_k$ and $\mathbf{D}_k$ store respectively the coefficients $A_{ik}, B_{ik}, C_{ik}$ and $D_{ik}$ with $k = 1, 2, \dots, N_h$ and $i = 1, 2, \dots, N$ are the Fourier coefficients to be found. Only the velocity and acceleration of the fluid-solid relative displacement is considered here because there is no stiffness term on $w$ in the equations of motion. Term-wise differentiation in time is achieved, such that

$$\dot{\mathbf{w}}_s(t) = \sum_{k=1}^{N_h} k\omega(-\mathbf{A}_k \sin k\omega t + \mathbf{B}_k \cos k\omega t), \quad \ddot{\mathbf{w}}_s(t) = \sum_{k=1}^{N_h} -k^2\omega^2(\mathbf{A}_k \cos k\omega t + \mathbf{B}_k \sin k\omega t), \quad (17a)$$

$$\ddot{\mathbf{w}}(t) = \sum_{k=1}^{N_h} k\omega(-\mathbf{C}_k \sin k\omega t + \mathbf{D}_k \cos k\omega t). \quad (17b)$$

The temporal contribution of the excitation force is assumed periodic and reads

$$\mathbf{q}_{\text{ex}}(t) = \frac{1}{2}\mathbf{A}_0^{\text{ex}} + \sum_{k=1}^{N_h} \mathbf{A}_k^{\text{ex}} \cos k\omega t + \mathbf{B}_k^{\text{ex}} \sin k\omega t, \quad (18)$$

where $\mathbf{A}_k^{\text{ex}}$ and $\mathbf{B}_k^{\text{ex}}$ are a priori known. Usually, the excitation force is considered monoharmonic and the only non-zero coefficient is $\mathbf{A}_1^{\text{ex}}$. The Fourier expansion of every entry of the nonlinear vector reads

$$f_i^{\text{NL}}(\dot{\mathbf{w}}(t)) = h\sigma_0\epsilon \left( \frac{1}{2}\Lambda_i + \sum_{k=1}^{N_h} (\Gamma_{ik} \cos k\omega t + \Pi_{ik} \sin k\omega t) \right), \quad (19)$$

where $\Lambda_i, \Gamma_{ik}$ and $\Pi_{ik}$ are nonlinear functions of the Fourier coefficients stored in $\mathbf{C}_k$ and $\mathbf{D}_k$, as follows:

$$\Lambda_i = \frac{1}{2} \int_S \left( \sum_{j=1}^{N} \Psi_j C_{j0} \right) \left| \sum_{r=1}^{N} \Psi_r C_{r0} \right| \Psi_i \, dx dy, \quad (20a)$$

$$\Gamma_{ik} = \frac{\omega}{\pi} \int_{0}^{2\pi/\omega} \int_{S} \left( \sum_{j=1}^{N} \Psi_j \sum_{\ell=1}^{N_h} (C_{j\ell} \cos \ell\omega\tau + D_{j\ell} \sin \ell\omega\tau) \right) \times \left| \sum_{r=1}^{N} \Psi_r \sum_{s=1}^{N_h} (C_{rs} \cos s\omega\tau + D_{rs} \sin s\omega\tau) \right| \cos k\omega\tau \Psi_i \, dx dy d\tau, \quad (20b)$$

$$\Pi_{ik} = \frac{\omega}{\pi} \int_{0}^{2\pi/\omega} \int_{S} \left( \sum_{j=1}^{N} \Psi_j \sum_{\ell=1}^{N_h} (C_{j\ell} \cos \ell\omega\tau + D_{j\ell} \sin \ell\omega\tau) \right) \times \left| \sum_{r=1}^{N} \Psi_r \sum_{s=1}^{N_h} (C_{rs} \cos s\omega\tau + D_{rs} \sin s\omega\tau) \right| \sin k\omega\tau \Psi_i \, dx dy d\tau. \quad (20c)$$

Substituting Equations (16), (17), (19) and (20) into Equation (10), and equating the coefficients associated with each harmonic component $\cos k\omega t$ and $\sin k\omega t$ with $k = 1, 2, \dots, N_h$, yields a system of $(4N) \times (2N_h+2)$ equations. In order to simplify the expressions, it is possible to rewrite Equation (10) in a matrix form in the frequency domain as

$$\mathcal{H}(\mathbf{x}, \omega) = \mathbf{G}(\omega)\mathbf{x}(\omega) + \tilde{\mathbf{f}}_{\text{NL}}(\mathbf{x}, \omega) - \tilde{\mathbf{f}}_{\text{ex}}(\omega) \quad (21)$$

where

$$\mathbf{x}(\omega) = \left[ \frac{1}{2}\mathbf{A}_0^\top, \frac{1}{2}\mathbf{C}_0^\top, \mathbf{A}_1^\top, \mathbf{B}_1^\top, \mathbf{C}_1^\top, \mathbf{D}_1^\top, \mathbf{A}_2^\top, \mathbf{B}_2^\top, \mathbf{C}_2^\top, \mathbf{D}_2^\top, \dots, \mathbf{A}_{N_h}^\top, \mathbf{B}_{N_h}^\top, \mathbf{C}_{N_h}^\top, \mathbf{D}_{N_h}^\top \right]^\top. \quad (22)$$
Based on the modal expansion and spatial semi-discretization procedure outlined above, Equation (10) can be expressed in matrix form:

$$ \mathbf{M}_1 \ddot{\mathbf{w}}_s(t) + \mathbf{M} \ddot{\mathbf{w}}(t) + \mathbf{K}_1 \mathbf{w}_s(t) = \mathbf{f}^{\text{ex}}(t), \quad (13a) $$

$$ \mathbf{M} \ddot{\mathbf{w}}_s(t) + \mathbf{M}_2 \ddot{\mathbf{w}}(t) + \mathbf{C} \dot{\mathbf{w}}(t) + \mathbf{f}^{\text{NL}}(\dot{\mathbf{w}}(t)) + \mathbf{K}_2 \mathbf{w}_s(t) = \mathbf{0}, \quad (13b) $$

where

$$ \begin{aligned} \mathbf{K}_1 &= D \int_S \nabla^2 \mathbf{\Psi}(x, y) \nabla^2 \mathbf{\Psi}^\top(x, y) \, dx \, dy, & \mathbf{K}_2 &= h \alpha M_f \int_S \nabla^2 \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) \, dx \, dy, \\ \mathbf{C} &= h \sigma_0 \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) \, dx \, dy, & \mathbf{M} &= h \rho_f \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) \, dx \, dy, \quad (14) \\ \mathbf{M}_1 &= h \rho_f \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) \, dx \, dy, & \mathbf{M}_2 &= \frac{h \rho_f \alpha_\infty}{\phi} \int_S \mathbf{\Psi}(x, y) \mathbf{\Psi}^\top(x, y) \, dx \, dy. \end{aligned} $$

where $S = L_x \times L_y$. Entry $i$ of the nonlinear force vector $\mathbf{f}^{\text{NL}}(\dot{\mathbf{w}}(t))$ reads

$$ f_i^{\text{NL}}(\dot{\mathbf{w}}(t)) = h \sigma_0 \xi \int_S \sum_{j=1}^{N} \dot{w}_j(t) \Psi_j(x, y) \left| \sum_{p=1}^{N} \dot{w}_p(t) \Psi_p(x, y) \right| \Psi_i(x, y) \, dx \, dy. \quad (15) $$

The Harmonic Balance Method (HBM) is implemented and the contributions are assumed to be periodic and therefore expressed as Fourier series truncated at the rank $N_h$ [34–36]:

$$ \mathbf{w}_s(t) = \frac{1}{2} \mathbf{A}_0 + \sum_{k=1}^{N_h} \mathbf{A}_k \cos k\omega t + \mathbf{B}_k \sin k\omega t \quad \text{and} \quad \dot{\mathbf{w}}(t) = \frac{1}{2} \mathbf{C}_0 + \sum_{k=1}^{N_h} \mathbf{C}_k \cos k\omega t + \mathbf{D}_k \sin k\omega t. \quad (16) $$

The vectors $\mathbf{A}_k$, $\mathbf{B}_k$, $\mathbf{C}_k$ and $\mathbf{D}_k$ store respectively the coefficients $A_{ik}$, $B_{ik}$, $C_{ik}$ and $D_{ik}$ with $k = 1, 2, \dots, N_h$ and $i = 1, 2, \dots, N$ are the Fourier coefficients to be found. Only the velocity and acceleration of the fluid-solid relative displacement is considered here because there is no stiffness term on $w$ in the equations of motion. Term-wise differentiation in time is achieved, such that

$$ \begin{aligned} \dot{\mathbf{w}}_s(t) &= \sum_{k=1}^{N_h} k\omega(-\mathbf{A}_k \sin k\omega t + \mathbf{B}_k \cos k\omega t), & \ddot{\mathbf{w}}_s(t) &= \sum_{k=1}^{N_h} -k^2\omega^2(\mathbf{A}_k \cos k\omega t + \mathbf{B}_k \sin k\omega t), \quad (17a) \\ \ddot{\mathbf{w}}(t) &= \sum_{k=1}^{N_h} k\omega(-\mathbf{C}_k \sin k\omega t + \mathbf{D}_k \cos k\omega t). & & \quad (17b) \end{aligned} $$

The temporal contribution of the excitation force is assumed periodic and reads

$$ \mathbf{q}^{\text{ex}}(t) = \frac{1}{2} \mathbf{A}_0^{\text{ex}} + \sum_{k=1}^{N_h} \mathbf{A}_k^{\text{ex}} \cos k\omega t + \mathbf{B}_k^{\text{ex}} \sin k\omega t, \quad (18) $$

where $\mathbf{A}_k^{\text{ex}}$ and $\mathbf{B}_k^{\text{ex}}$ are a priori known. Usually, the excitation force is considered monoharmonic and the only non-zero coefficient is $\mathbf{A}_1^{\text{ex}}$. The Fourier expansion of every entry of the nonlinear vector reads

$$ f_i^{\text{NL}}(\dot{\mathbf{w}}(t)) = h \sigma_0 \xi \left( \frac{1}{2} \Lambda_i + \sum_{k=1}^{N_h} (\Gamma_{ik} \cos k\omega t + \Pi_{ik} \sin k\omega t) \right), \quad (19) $$

where $\Lambda_i$, $\Gamma_{ik}$ and $\Pi_{ik}$ are nonlinear functions of the Fourier coefficients stored in $\mathbf{C}_k$ and $\mathbf{D}_k$, as follows:

$$ \begin{aligned} \Lambda_i &= \frac{1}{2} \int_S \left( \sum_{j=1}^N \Psi_j C_{j0} \right) \left| \sum_{r=1}^N \Psi_r C_{r0} \right| \Psi_i \, dx \, dy, & (20a) \\ \Gamma_{ik} &= \frac{\omega}{\pi} \int_0^{2\pi/\omega} \int_S \left( \sum_{j=1}^N \Psi_j \sum_{\ell=1}^{N_h} (C_{j\ell} \cos \ell\omega\tau + D_{j\ell} \sin \ell\omega\tau) \right) \times \left| \sum_{r=1}^N \Psi_r \sum_{s=1}^{N_h} (C_{rs} \cos s\omega\tau + D_{rs} \sin s\omega\tau) \right| \cos k\omega\tau \Psi_i \, dx \, dy \, d\tau, & (20b) \end{aligned} $$

$$ \begin{aligned} \Pi_{ik} &= \frac{\omega}{\pi} \int_0^{2\pi/\omega} \int_S \left( \sum_{j=1}^N \Psi_j \sum_{\ell=1}^{N_h} (C_{j\ell} \cos \ell\omega\tau + D_{j\ell} \sin \ell\omega\tau) \right) \times \left| \sum_{r=1}^N \Psi_r \sum_{s=1}^{N_h} (C_{rs} \cos s\omega\tau + D_{rs} \sin s\omega\tau) \right| \sin k\omega\tau \Psi_i \, dx \, dy \, d\tau. & (20c) \end{aligned} $$
stores the Fourier's coefficients $A_{ik}$, $B_{ik}$, $C_{ik}$ and $D_{ik}$ (the $\omega$ contributions are dropped for readability purpose). The matrix $G = \text{blockdiag}(0, G_1, G_2, \dots, G_{N_h})$ is a block matrix with each block of dimension of $4N$. The elementary block $k$ reads

$$G_k(\omega) = \begin{bmatrix} -\omega^2 k^2 \mathbf{M}_1 + \mathbf{K}_1 & \mathbf{0} & \mathbf{0} & \omega k \mathbf{M} \\ \mathbf{0} & -\omega^2 k^2 \mathbf{M}_1 + \mathbf{K}_1 & -\omega k \mathbf{M} & \mathbf{0} \\ -\omega^2 k^2 \mathbf{M}_2 + \mathbf{K}_2 & \mathbf{0} & \mathbf{C} & \omega k \mathbf{M} \\ \mathbf{0} & -\omega^2 k^2 \mathbf{M}_2 + \mathbf{K}_2 & -\omega k \mathbf{M} & \mathbf{C} \end{bmatrix}. \quad (23)$$

The $n(4N_h + 2)$-dimensional vectors $\tilde{\mathbf{f}}_{\text{ex}}(\omega)$ and $\tilde{\mathbf{f}}_{\text{NL}}(\mathbf{x}, \omega)$ contain respectively the Fourier coefficients of the external forcing contribution and those of the nonlinear force contribution, and read

$$\tilde{\mathbf{f}}_{\text{ex}}(\omega) = \left[ \frac{1}{2} (\mathbf{A}_0^{\text{ex}})^{\top}, \mathbf{0}, (\mathbf{A}_1^{\text{ex}})^{\top}, (\mathbf{B}_1^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{A}_2^{\text{ex}})^{\top}, (\mathbf{B}_2^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0}, \dots, (\mathbf{A}_{N_h}^{\text{ex}})^{\top}, (\mathbf{B}_{N_h}^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0} \right]^{\top} \quad (24a)$$

$$\tilde{\mathbf{f}}_{\text{NL}}(\mathbf{x}, \omega) = h\sigma_0\epsilon \left[ \mathbf{0}, \mathbf{\Lambda}^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_1)^{\top}, (\mathbf{\Pi}_1)^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_2)^{\top}, (\mathbf{\Pi}_2)^{\top}, \dots, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_{N_h})^{\top}, (\mathbf{\Pi}_{N_h})^{\top} \right]^{\top}. \quad (24b)$$

The vectors $\mathbf{\Lambda}$, $\mathbf{\Gamma}_k(\omega)$ and $\mathbf{\Pi}_k(\omega)$ are obtained for each value of $\omega$ through Equation (20). For known values of $\mathbf{A}_k^{\text{ex}}$ and $\mathbf{B}_k^{\text{ex}}$, Equation (21) is solved using the Python built-in `fsolve` function, which returns the roots of the nonlinear equation $\mathcal{H}(\mathbf{x}, \omega) = 0$ and leads to the unknown coefficients $\mathbf{A}_k$, $\mathbf{B}_k$, $\mathbf{C}_k$ and $\mathbf{D}_k$. The `fsolve` function is a version of the root-solving Powell's hybrid method. This approach combines the Gauss-Newton algorithm with gradient descent and approximates the Jacobian matrix by direct differences at the starting point. This method causes no convergence issue even with the nonlinear absolute value function involved in the formulation of the problem, meaning there was no reason to opt for a potentially more advanced solver.

Using the single-mode assumption, the temporal and spatial contributions in Equation (15) can be decoupled. Otherwise, the quadratic damping term implies elastic coupling between the modes due to the absolute value function in Equation (20). Then, temporal and spatial dependencies cannot be separated and Equation (20) must be calculated for each nonlinear step of the algorithm. The integral is approximated numerically using the Python built-in function `nquad`.

## 3.2. Convergence analysis

First, a convergence analysis is performed on the number of harmonics and the number of dof in the approximation. The overall residual is denoted $r(\omega) = \|\mathbf{r}(\omega)\|$ with $\mathbf{r}(\omega) = [(\mathbf{r}^{(1)}(\omega))^{\top}, (\mathbf{r}^{(2)}(\omega))^{\top}]^{\top}$ where the vectors $\mathbf{r}^{(1)}(\omega)$ and $\mathbf{r}^{(2)}(\omega)$ are determined from Equations (13a) and (13b) by the integrals

$$\mathbf{r}^{(1)}(\omega) = \frac{\omega}{2\pi} \int_{0}^{2\pi/\omega} \left( \mathbf{M}_1 \ddot{\mathbf{w}}_s(t) + \mathbf{M} \ddot{\mathbf{w}}(t) + \mathbf{K}_1 \mathbf{w}_s(t) - \mathbf{f}_{\text{ex}}(t) \right) dt, \quad (25a)$$

$$\mathbf{r}^{(2)}(\omega) = \frac{\omega}{2\pi} \int_{0}^{2\pi/\omega} \left( \mathbf{M} \ddot{\mathbf{w}}_s(t) + \mathbf{M}_2 \ddot{\mathbf{w}}(t) + \mathbf{C} \ddot{\mathbf{w}}(t) + \mathbf{f}_{\text{NL}}(\dot{\mathbf{w}}(t)) + \mathbf{K}_2 \mathbf{w}_s(t) \right) dt. \quad (25b)$$

Equation (25a) is numerically computed using the Python built-in `nquad`. The vector $\mathbf{f}_{\text{NL}}(\dot{\mathbf{w}}(t))$ depends on the mode shape and therefore on $x$ and $y$. After discretizing the three axes $x$, $y$ and $t$, a three-dimensional Riemann sum is used to numerically approximate Equation (25b) thanks to a Gauss-Legendre quadrature with 100 points for spatial discretization and 300 points for temporal discretization to guarantee accuracy. The vectors $\mathbf{w}_s(t)$ and $\mathbf{w}(t)$ are obtained by injecting the Fourier coefficients stored in $\mathbf{x}(\omega)$ in Equations (17), (19) and (20). The residual and its corresponding computation time $t_f$ are evaluated for a simply-supported aluminum MPP of dimension 490 mm × 570 mm × 1 mm with $d = 2.2$ mm, $\phi = 10\%$ and $\epsilon = 0.87$ s m$^{-1}$. The sensitivity to $N$ and $N_h$ of the normalized residual

$$L_r(\omega) = 20 \log_{10} \left( \frac{r(\omega)}{r_{\text{ref}}(\omega)} \right), \quad (26)$$

is indicated in Figure 1 for $\omega = \omega_1$ where $\omega_1$ is the first resonance frequency of the linear MPP. The solution vector $\mathbf{x}(\omega)$ is obtained for $\omega_1 = 115$ rad s$^{-1}$ from the procedure presented in Section 3. The residual is computed for $N = 1, 3, 5, 7$ and for $N_h = 1, 3, 5, 9, 13, 17$. Even numbers of harmonics are not studied, as the nonlinearity is odd, so only odd harmonics contribute. The residual reference $r_{\text{ref}}(\omega_1)$ in Equation (26) is set to $1 \times 10^{-11}$, which corresponds to $r(\omega_1)$ obtained for $N = 7$ and $N_h = 17$ and to the relative error between two consecutive iterations imposed as a stopping criterion in `fsolve`. In light of Figure 1, it is clear that for moderate excitations, $N_h = 1$ is sufficient to accurately predict the structural response of the system. The residual for $N = N_h = 1$ is equal to $r(\omega_1) = 4 \times 10^{-7}$. Moreover, the computation time increases significantly as $N$ and $N_h$ increase, from 0.64 s for $N = N_h = 1$ to around 3 hours for $N = 7$ and $N_h = 17$. Most of the computing time is allocated to calculating triple integrals (Equation (20)). In fact, the time required to compute an integral of Equation (20) increases from 0.02 s for $N = N_h = 1$ to 977 s for $N = 7$ and $N_h = 17$. Consequently, in the remainder of this paper, only the fundamental harmonic is considered in the HBM resolution.
Substituting Equations (16), (17), (19) and (20) into Equation (10), and equating the coefficients associated with each harmonic component cos $k\omega t$ and sin $k\omega t$ with $k = 1, 2, \dots, N_h$, yields a system of $(4N) \times (2N_h+2)$ equations. In order to simplify the expressions, it is possible to rewrite Equation (10) in a matrix form in the frequency domain as

$$ \mathcal{H}(\mathbf{x}, \omega) = \mathbf{G}(\omega)\mathbf{x}(\omega) + \tilde{\mathbf{f}}^{\text{NL}}(\mathbf{x}, \omega) - \tilde{\mathbf{f}}^{\text{ex}}(\omega) \quad (21) $$

where

$$ \mathbf{x}(\omega) = \left[ \frac{1}{2} \mathbf{A}_0^\top, \frac{1}{2} \mathbf{C}_0^\top, \mathbf{A}_1^\top, \mathbf{B}_1^\top, \mathbf{C}_1^\top, \mathbf{D}_1^\top, \mathbf{A}_2^\top, \mathbf{B}_2^\top, \mathbf{C}_2^\top, \mathbf{D}_2^\top, \dots, \mathbf{A}_{N_h}^\top, \mathbf{B}_{N_h}^\top, \mathbf{C}_{N_h}^\top, \mathbf{D}_{N_h}^\top \right]^\top. \quad (22) $$

stores the Fourier's coefficients $A_{ik}$, $B_{ik}$, $C_{ik}$ and $D_{ik}$ (the $\omega$ contributions are dropped for readability purpose). The matrix $\mathbf{G} = \text{blockdiag}(\mathbf{0}, \mathbf{G}_1, \mathbf{G}_2, \dots, \mathbf{G}_{N_h})$ is a block matrix with each block of dimension of $4N$. The elementary block $k$ reads

$$ \mathbf{G}_k(\omega) = \begin{bmatrix} -\omega^2 k^2 \mathbf{M}_1 + \mathbf{K}_1 & \mathbf{0} & \mathbf{0} & \omega k \mathbf{M} \\ \mathbf{0} & -\omega^2 k^2 \mathbf{M}_1 + \mathbf{K}_1 & -\omega k \mathbf{M} & \mathbf{0} \\ -\omega^2 k^2 \mathbf{M}_2 + \mathbf{K}_2 & \mathbf{0} & \mathbf{C} & \omega k \mathbf{M} \\ \mathbf{0} & -\omega^2 k^2 \mathbf{M}_2 + \mathbf{K}_2 & -\omega k \mathbf{M} & \mathbf{C} \end{bmatrix}. \quad (23) $$

The $n(4N_h + 2)$-dimensional vectors $\tilde{\mathbf{f}}^{\text{ex}}(\omega)$ and $\tilde{\mathbf{f}}^{\text{NL}}(\mathbf{x}, \omega)$ contain respectively the Fourier coefficients of the external forcing contribution and those of the nonlinear force contribution, and read

$$ \tilde{\mathbf{f}}^{\text{ex}}(\omega) = \left[ \frac{1}{2}(\mathbf{A}_0^{\text{ex}})^{\top}, \mathbf{0}, (\mathbf{A}_1^{\text{ex}})^{\top}, (\mathbf{B}_1^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{A}_2^{\text{ex}})^{\top}, (\mathbf{B}_2^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0}, \dots, (\mathbf{A}_{N_h}^{\text{ex}})^{\top}, (\mathbf{B}_{N_h}^{\text{ex}})^{\top}, \mathbf{0}, \mathbf{0} \right]^{\top} \quad (24a) $$

$$ \tilde{\mathbf{f}}^{\text{NL}}(\mathbf{x}, \omega) = h\sigma_0\xi [\mathbf{0}, \mathbf{\Lambda}^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_1)^{\top}, (\mathbf{\Pi}_1)^{\top}, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_2)^{\top}, (\mathbf{\Pi}_2)^{\top}, \dots, \mathbf{0}, \mathbf{0}, (\mathbf{\Gamma}_{N_h})^{\top}, (\mathbf{\Pi}_{N_h})^{\top}]^{\top}. \quad (24b) $$

The vectors $\mathbf{A}$, $\mathbf{\Gamma}_k(\omega)$ and $\mathbf{\Pi}_k(\omega)$ are obtained for each value of $\omega$ through Equation (20). For known values of $\mathbf{A}_k^{\text{ex}}$ and $\mathbf{B}_k^{\text{ex}}$, Equation (21) is solved using the Python built-in `fsolve` function, which returns the roots of the nonlinear equation $\mathcal{H}(\mathbf{x}, \omega) = 0$ and leads to the unknown coefficients $\mathbf{A}_k$, $\mathbf{B}_k$, $\mathbf{C}_k$ and $\mathbf{D}_k$. The `fsolve` function is a version of the root-solving Powell's hybrid method. This approach combines the Gauss-Newton algorithm with gradient descent and approximates the Jacobian matrix by direct differences at the starting point. This method causes no convergence issue even with the nonlinear absolute value function involved in the formulation of the problem, meaning there was no reason to opt for a potentially more advanced solver.

Using the single-mode assumption, the temporal and spatial contributions in Equation (15) can be decoupled. Otherwise, the quadratic damping term implies elastic coupling between the modes due to the absolute value function in Equation (20). Then, temporal and spatial dependencies cannot be separated and Equation (20) must be calculated for each nonlinear step of the algorithm. The integral is approximated numerically using the Python built-in function `nquad`.

## 3.2. Convergence analysis

First, a convergence analysis is performed on the number of harmonics and the number of dof in the approximation. The overall residual is denoted $r(\omega) = \|\mathbf{r}(\omega)\|$ with $\mathbf{r}(\omega) = [(\mathbf{r}^{(1)}(\omega))^{\top}, (\mathbf{r}^{(2)}(\omega))^{\top}]^{\top}$ where the vectors $\mathbf{r}^{(1)}(\omega)$ and $\mathbf{r}^{(2)}(\omega)$ are determined from Equations (13a) and (13b) by the integrals

$$ \mathbf{r}^{(1)}(\omega) = \frac{\omega}{2\pi} \int_{0}^{2\pi/\omega} \left( \mathbf{M}_1 \ddot{\mathbf{w}}_s(t) + \mathbf{M} \ddot{\mathbf{w}}(t) + \mathbf{K}_1 \mathbf{w}_s(t) - \mathbf{f}^{\text{ex}}(t) \right) dt, \quad (25a) $$

$$ \mathbf{r}^{(2)}(\omega) = \frac{\omega}{2\pi} \int_{0}^{2\pi/\omega} \left( \mathbf{M} \ddot{\mathbf{w}}_s(t) + \mathbf{M}_2 \ddot{\mathbf{w}}(t) + \mathbf{C} \ddot{\mathbf{w}}(t) + \mathbf{f}^{\text{NL}}(\dot{\mathbf{w}}(t)) + \mathbf{K}_2 \mathbf{w}_s(t) \right) dt. \quad (25b) $$

Equation (25a) is numerically computed using the Python built-in `nquad`. The vector $\mathbf{f}^{\text{NL}}(\dot{\mathbf{w}}(t))$ depends on the mode shape and therefore on $x$ and $y$. After discretizing the three axes $x$, $y$ and $t$, a three-dimensional Riemann sum is used to numerically approximate Equation (25b) thanks to a Gauss-Legendre quadrature with 100 points for spatial discretization and 300 points for temporal discretization to guarantee accuracy. The vectors $\mathbf{w}_s(t)$ and $\mathbf{w}(t)$ are obtained by injecting the Fourier coefficients stored in $\mathbf{x}(\omega)$ in Equations (17), (19) and (20). The residual and its corresponding computation time $t_f$ are evaluated for a simply-supported aluminum MPP of dimensions 490 mm × 570 mm × 1 mm with $d = 2.2$ mm, $\phi = 10\%$ and $\xi = 0.87$ s m$^{-1}$. The sensitivity to $N$ and $N_h$ of the normalized residual

$$ L_r(\omega) = 20 \log_{10} \left( \frac{r(\omega)}{r_{\text{ref}}(\omega)} \right), \quad (26) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Influence of the number of harmonics on the forced response of the first mode of a simply supported microperforated plate. (a) — normalized residual calculated from Equation (25) for $F_{\text{ext}} = 0.1$ N; (b) — computation time (inverted axes for readability).</figcaption></figure>

# 4. Sensitivity of microperforated plates vibratory response to nonlinear damping

## 4.1. Frequency response

Equation (13) is solved using the methodology described in Section 3 considering a single mode and only the contribution of the fundamental harmonic that is $N = N_h = 1$. In Figure 2, the mobility at the center of an MPP is plotted to investigate the influence of the amplitude of the external excitation on the vibration response. An aluminum simply supported microperforated plate is excited by a point force at $(x_e, y_e) = (60 \text{ mm}, 70 \text{ mm})$. The mobility at point $(x_g, y_g) = (L_x/2, L_y/2)$, that is the solid displacement normalized by the amplitude of the external force, is plotted in function of the normalized angular pulsation of the excitation $\omega/\omega_1$. Here $\omega_1$ is the angular frequency of the first linear resonance of the MPP with added damping, obtained by assuming $\varepsilon = 0$ in Equation (10). The mobility is plotted for different amplitudes of $F_{\text{ext}} \in [10^{-4}; 0.4]$ N. The aluminum microperforated plate is simply supported with $L_x \times L_y \times h = 490 \text{ mm} \times 570 \text{ mm} \times 1 \text{ mm}$, $d = 2.2 \text{ mm}$ and $\phi = 10\%$. The corresponding Forchheimer's parameter is $\varepsilon = 0.87 \text{ s m}^{-1}$. In this configuration, the characteristic frequency is lower than the first MPP resonance frequency. Setting $\varepsilon$ to zero and solving Equations (10a) and (10b) leads to the MPP

<figure><img src="image_4.png" /><figcaption>Figure 2: Mobility at the center of the simply-supported MPP for different amplitudes of the external force with $d = 2.2 \text{ mm}$ and $\phi = 10\%$. The mobility is given for $F_{\text{ext}} \in [10^{-4}; 0.4]$ N and for: (- - -) $F_{\text{ext}} = 10^{-4}$ N that correspond to the linear case; (—•—) $F_{\text{ext}} = 0.16$ N which corresponds to the case where the maximum of the normalized displacement passes through a minimum. The $\varepsilon$-parameter is obtained from Equation (9). The resolution procedure is performed for $N = N_h = 1$.</figcaption></figure>

response in the linear case, and thus to $\omega_1$.

In Figure 2, the linear case is obtained for $F_{\text{ext}} = 10^{-4}$ N. In this configuration, the first resonance frequency is $\omega_1$ and the displacement amplitude is governed by the added viscous damping in the linear regime and thus by the loss factor given in Equation (4). When the excitation amplitude $F_{\text{ext}}$ increases, there are two prominent effects: the equivalent stiffness of the system is reduced, and the added damping can increase or decrease with the excitation amplitude.
is indicated in Figure 1 for $\omega = \omega_1$ where $\omega_1$ is the first resonance frequency of the linear MPP. The solution vector $\mathbf{x}(\omega)$ is obtained for $\omega_1 = 115 \text{ rad s}^{-1}$ from the procedure presented in Section 3. The residual is computed for $N = 1, 3, 5, 7$ and for $N_h = 1, 3, 5, 9, 13, 17$. Even numbers of harmonics are not studied, as the nonlinearity is odd, so only odd harmonics contribute. The residual reference $r_{\text{ref}}(\omega_1)$ in Equation (26) is set to $1 \times 10^{-11}$, which corresponds to $r(\omega_1)$ obtained for $N = 7$ and $N_h = 17$ and to the relative error between two consecutive iterations imposed as a stopping criterion in `fsolve`. In

<figure><img src="image_3.png" /><figcaption>Figure 1: Influence of the number of harmonics on the forced response of the first mode of a simply supported microperforated plate. (a) normalized residual calculated from Equation (25) for $F_{\text{ext}} = 0.1 \text{ N}$; (b) computation time (inverted axes for readability).</figcaption></figure>

light of Figure 1, it is clear that for moderate excitation, $N_h = 1$ is sufficient to accurately predict the structural response of the system. The residual for $N = N_h = 1$ is equal to $r(\omega_1) = 4 \times 10^{-7}$. Moreover, the computation time increases significantly as $N$ and $N_h$ increase, from 0.64 s for $N = N_h = 1$ to around 3 hours for $N = 7$ and $N_h = 17$. Most of the computing time is allocated to calculating triple integrals (Equation (20)). In fact, the time required to compute an integral of Equation (20) increases from 0.02 s for $N = N_h = 1$ to 977 s for $N = 7$ and $N_h = 17$. Consequently, in the remainder of this paper, only the fundamental harmonic is considered in the HBM resolution.

# 4. Sensitivity of microperforated plates vibratory response to nonlinear damping

## 4.1. Frequency response

Equation (13) is solved using the methodology described in Section 3 considering a single mode and only the contribution of the fundamental harmonic that is $N = N_h = 1$. In Figure 2, the mobility at the center of an MPP is plotted to investigate the influence of the amplitude of the external excitation on the vibration response. An aluminum simply supported microperforated plate is excited by a point force at $(x_e, y_e) = (60 \text{ mm}, 70 \text{ mm})$. The mobility at point $(x_g, y_g) = (L_x/2, L_y/2)$, that is the solid displacement normalized by the amplitude of the external force, is plotted in function of the normalized angular pulsation of the excitation $\omega/\omega_1$. Here $\omega_1$ is the angular frequency of the first linear resonance of the MPP with added damping, obtained by assuming $\xi = 0$ in Equation (10). The mobility is plotted for different amplitudes of $F_{\text{ext}} \in [10^{-4}; 1] \text{ N}$. The aluminum microperforated plate is simply supported with $L_x \times L_y \times h = 490 \text{ mm} \times 570 \text{ mm} \times 1 \text{ mm}$, $d = 2.2 \text{ mm}$ and $\phi = 10 \%$. The corresponding Forchheimer's parameter is $\xi = 0.87 \text{ s m}^{-1}$. In this configuration, the characteristic frequency is lower than the first MPP resonance frequency. Setting $\xi$ to zero and solving Equations (10a) and (10b) leads to the MPP response in the linear case, and thus to $\omega_1$.

In Figure 2, the linear case is obtained for $F_{\text{ext}} = 10^{-4} \text{ N}$. In this configuration, the first resonance frequency is $\omega_1$ and the displacement amplitude is governed by the added viscous damping in the linear regime and thus by the loss factor given in Equation (4). When the excitation amplitude $F_{\text{ext}}$ increases, there are two prominent effects: the equivalent stiffness of the system is reduced, and the added damping can increase or decrease with the excitation amplitude.

In fact, the maximum displacement amplitude decreases and passes through a minimum, which corresponds to $F_{\text{ext}} \approx 0.2 \text{ N}$ in this particular case, before rising with $F_{\text{ext}}$. This minimum is reached for a specific value of the driving frequency, referred to as the characteristic nonlinear frequency $\omega_c^{\text{NL}}$ and discussed further in the next section. As mentioned in Section 2.2, the linear and, by extension, nonlinear damping added by microperforations is related to the thickness of the boundary layers. For a microperforated plate, viscous dissipation and thus the damping is governed by the airflow resistivity. Due to the Forchheimer correction presented in Equation (8), $\sigma$ depends on the relative fluid-solid velocity. Thus, viscous dissipation can reach a maximum at a particular relative fluid-solid velocity, referred to in
In fact, the maximum displacement amplitude decreases and passes through a minimum, which corresponds to $F_{ext} \approx 0.16$ N in this particular case, before rising with $F_{ext}$. As mentioned in Section 2.2, the linear and, by extension, nonlinear damping added by microperforations is related to the thickness of the boundary layers. For a microperforated plate, viscous dissipation and thus the damping is governed by the airflow resistivity. Due to the Forchheimer correction presented in Equation (8), $\sigma$ depends on the relative fluid-solid velocity. Thus, viscous dissipation can reach a maximum at a particular relative fluid-solid velocity, referred to in this paper as critical relative velocity in reference to the critical mach number defined by [7]. This critical velocity is a function of the perforation parameters<fnref n="3" />. As the relative fluid-solid velocity is linked to the excitation amplitude $F_{ext}$, there is an excitation amplitude value for which the size of the boundary layers maximizes the viscous damping effect. In the zero forcing limit, the plate velocity and the relative fluid-solid velocity are very low. Resistive effects are equal to $\sigma_0$ and do not depend on external excitation. This corresponds to the linear case represented by the dashed line in Figure 2. In the limit of infinite forcing, solid and fluid velocities will also tend towards infinity. In this case, the resistivity, which varies linearly with the relative fluid-solid velocity, tends toward infinity, and no dissipative effects occur in the perforations. In this case, the added damping will not be maximal and the amplitude of the peaks will increase with forcing. Between these two limits, damping passes through a maximum<fnref n="4" />.

## 4.2. Normalized resistance and maximum of added damping in nonlinear regime

It has been established that in a linear context, the maximum added damping is reached for $f = f_c$ [5], but this maximum evolves with relative velocity under nonlinear regime. The evolution of the maximum of the added damping with relative velocity is analyzed using the normalized resistance. In Figure 2, it was observed that the added damping reaches a maximum for $F_{ext} = 0.16$ N (this corresponds to the minimum of the amplitude maximum in Figure 2). In this configuration, the amplitude of the relative fluid-solid velocity is equal to $2.4$ m s⁻¹ and the normalized nonlinear resistance [18]

$$R_{NL}(\dot{w}) = R_L(1 + \varepsilon|\dot{w}|), \qquad (27)$$

with the normalized linear resistance

$$R_L = \frac{h\sigma_0}{\rho_f c_0 \phi^2} \approx 1. \qquad (28)$$

In Equation (28), $c_0$ is the sound velocity in the fluid. The linear resistance depends on fluid parameters ($c_0$ and $\rho_f$) and perforation parameters ($\phi$, $d$ and $h$).

In an acoustics framework, it was noted that the maximum of the acoustic absorption coefficient could pass through a maximum for an MPP subjected to a high level of fluid pressure. From the results of Tayong et al. [7], it can be observed that the maximum of the acoustic absorption coefficient for an MPP, backed by an air cavity, is obtained when the imaginary part of the total normalized impedance is equal to zero and the real part is equal to 1⁵. For an MPP, the dissipation mechanisms are the same for acoustic absorption and vibration damping. In both cases, viscous effects are maximized when normalized resistance is equal to 1. In the acoustic nonlinearity regime, the maximum of added damping can only be reached if $R_{NL} = 1$. According to Equation (27), the nonlinear normalized resistance is an increasing function of $|\dot{w}|$. It can only increase with the relative fluid-solid velocity. Depending on the value of $R_L$, three cases can be distinguished:

1. $R_L < 1$ — The maximum viscous dissipation is not reached in the linear regime but can be reached in the nonlinear regime. Indeed, the maximum of the added damping can pass through a maximum by increasing the amplitude of the external forcing $F_{ext}$ and the corresponding magnitude of $\dot{w}$. Thus, $R_{NL} = 1$ is obtained for a particular value of $|\dot{w}|$ called critical relative velocity. Consequently, the maximum of added damping obtained in a linear case for $f = f_c$ increase until $V_{cr}$ which correspond to its maximum and then decrease according to the relative fluid-solid velocity.

2. $R_L = 1$ — This corresponds to the limit case. The maximum of the maximum added damping is reached in the linear regime. The maximum of added damping decreases with relative fluid-solid velocity from its maximum value. The dissipation mechanisms are maximal for both linear and nonlinear regimes.

<footnote n="3">³Details concerning this critical velocity are given in Section 4.3.</footnote>

<footnote n="4">⁴Similar behaviours have been observed for systems with dry friction where the linear term is proportional to $|\dot{u}|$ with $\dot{u}$ the structural velocity [37, 38]. However, for MPP, it is important to emphasize that these observations are only valid under a regime of nonlinearity exclusively linked to the fluid; all other nonlinearities (in particular those related to the solid) are neglected.</footnote>

<footnote n="5">⁵The reader is invited to refer to Equation (19) of [7] for further details.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: Mobility at the center of the simply-supported MPP for different amplitudes of the external force with $d = 2.2$ mm and $\phi = 10\%$. The mobility is given for $F_{\text{ext}} \in [10^{-4}; 1]$ N and for: (---) $F_{\text{ext}} = 10^{-4}$ N that corresponds to the linear case; (---) $F_{\text{ext}} = 0.2$ N which corresponds to the case where the maximum of the normalized displacement passes through a minimum. The $\xi$-parameter is obtained from Equation (9). The resolution procedure is performed for $N = N_h = 1$.</figcaption></figure>

this paper as critical velocity<fnref n="4" /> in reference to the critical Mach number defined by [7]. This critical velocity is a function of the perforation parameters. As the relative fluid-solid velocity is linked to the excitation amplitude $F_{\text{ext}}$, there is an excitation amplitude value for which the size of the boundary layers maximizes the viscous damping effect. In the zero forcing limit, the plate velocity and the relative fluid-solid velocity are very low. Resistive effects are equal to $\sigma_0$ and do not depend on external excitation. This corresponds to the linear case represented by the dashed line in Figure 2. In the limit of infinite forcing, solid and fluid velocities will also tend towards infinity. In this case, the resistivity, which varies linearly with the relative fluid-solid velocity, tends toward infinity, and no dissipative effects occur in the perforations. In this case, the added damping will not be maximal and the amplitude of the peaks will increase with forcing. Between these two limits, damping passes through a maximum. Note that similar behaviours have been observed for systems with dry friction where the linear term is proportional to $|\dot{u}|$ with $\dot{u}$ the structural velocity [37, 38] or for the system called Lanchester dampers [39, 40]. However, for MPP, it is important to emphasize that these observations are only valid under a regime of nonlinearity exclusively linked to the fluid; all other nonlinearities (in particular those related to the solid) are neglected.

## 4.2. Maximum of added damping in nonlinear regime

It has been established that, in a linear context, the maximum added damping is achieved when $\omega = \omega_c$ [5], where $\omega$ is the driving angular frequency. When the ratio $\omega/\omega_c$ differs from 1, the maximum added damping cannot be attained in the linear regime.

In the nonlinear regime, the added damping depends on the amplitude of the external excitation and, consequently, on the relative fluid-solid velocity (cf. Figure 2). To investigate the evolution of the added damping with respect to the relative fluid-solid velocity, we propose adapting Equation (4) to the acoustic nonlinear framework by modifying the airflow resistivity in accordance with the Forchheimer law described in Equation (8). Since $N = N_h = 1$, the time and space dependencies in $\dot{w}$ are separable, allowing it to be expressed in complex form as

$$ \dot{w}(x, y, t) = \Psi_1(x, y)V \exp(j\omega t), \quad V = \sqrt{C_{11}^2 + D_{11}^2}, \qquad (27) $$

where $C_{11}$ and $D_{11}$ are the unknown coefficients obtained from HBM. By substituting Equation (27) into Equation (8), the resistivity becomes a function of $V$, which is then inserted back into Equation (4). As a result, the loss factor for the nonlinear case is obtained, here only for the first mode and its fundamental harmonic, as:

$$ \eta_1(\omega, V) = -\frac{h\omega I_1^{(2)} I_1^{(3)}}{D I_1^{(1)}} \cdot \frac{\alpha M_f \phi^2 \rho_f \sigma_0 (I_1^{(3)} + \xi V I_1^{(4)})}{(\alpha_\infty \rho_f \omega I_1^{(3)})^2 + (\phi \sigma_0 (I_1^{(3)} + \xi V I_1^{(4)}))^2}, \qquad (28) $$

<footnote n="4">⁴Details concerning the critical velocity are given in Section 4.2.</footnote>
3. $R_L > 1$ — The normalized nonlinear resistance will always be greater than 1, and the maximum of the maximum added damping will not be reached in either linear or nonlinear regime. The plate parameters cannot achieve maximum viscous effects in the nonlinear regime.

In the following, the preceding paragraphs are discussed using an example. To this end, three MPPs whose perforation parameters have been chosen to illustrate the following three cases are analyzed: $R_L > 1$ for MPP①, $R_L = 1$ for MPP②, and $R_L < 1$ for MPP③. The perforation ratio for all MPP considered is set to $\phi = 10\%$. The linear resistance depending on $\sigma_0$ the perforation diameter is adjusted to accommodate the three required cases. For each MPP, the thickness is fixed to $h = 1$ mm and the length and width are adjusted so that the first resonance frequency is close to the characteristic frequency of each MPP. The MPP parameters, as well as $f_c$ and $R_L$ are stored in Table 1. For all MPP, the Forchheimer's

<table><caption>Table 1: Microperforated plate parameters for three cases: $R_L > 1$ for MPP①, $R_L = 1$ for MPP②, $R_L < 1$ for MPP③. For the three MPP, $L_x$ and $L_y$ are chosen so that the first resonant frequency of the structure coincides with $f_c$.</caption><thead><tr><th></th><th>MPP①</th><th>MPP②</th><th>MPP③</th></tr></thead><tbody><tr><td>$L_x$</td><td>0.248</td><td>0.405</td><td>0.800</td></tr><tr><td>$L_y$</td><td>0.300</td><td>0.406</td><td>0.697</td></tr><tr><td>$d$ [$\times 10^{-3}$]</td><td>0.800</td><td>1.195</td><td>2.20</td></tr><tr><td>$\varepsilon$</td><td>0.12</td><td>0.26</td><td>0.87</td></tr><tr><td>$f_c$</td><td>84.1</td><td>30.4</td><td>7.3</td></tr><tr><td>$R_L$</td><td>2.23</td><td>1</td><td>0.29</td></tr></tbody></table>

parameter is obtained from Equation (9).

In order to represent the added damping in the nonlinear regime, Equation (4) is adapted to the acoustic nonlinear framework by modifying the airflow resistivity using the Forchheimer law given in Equation (8). Again, since $N = N_h = 1$, the time and space dependencies in $\dot{w}$ are separable and it is possible to express it in a complex form as

$$ \dot{w}(x, y, t) = \Psi_1(x, y)V \exp(j\omega t), \quad V = \sqrt{C_{11}^2 + D_{11}^2}, \qquad (29) $$

where $C_{11}$ and $D_{11}$ are the unknown coefficients obtained from HBM. Substituting Equation (29) in Equation (4) yields an expression of the loss factor in the nonlinear case only for the first mode with the sole fundamental harmonic:

$$ \eta_1(\omega, \varepsilon V) = - \frac{h\omega I_1^{(2)} I_1^{(3)}}{D I_1^{(1)}} \cdot \frac{\alpha M_f \phi^2 \rho_f \sigma_0 (I_1^{(3)} + \varepsilon V I_1^{(4)})}{(\alpha_\infty \rho_f \omega I_1^{(3)})^2 + (\phi \sigma_0 (I_1^{(3)} + \varepsilon V I_1^{(4)}))^2}, \qquad (30) $$

where

$$ I_1^{(4)} = \int_S \Psi_1(x, y)^2 |\Psi_1(x, y)| dx dy \qquad (31) $$

carries the Forchheimer acoustic nonlinear mechanism. The added loss factor is plotted as a function of the amplitude of the relative fluid-solid velocity for MPP①, MPP② and MPP③ in Figure 3. The studied MPP are excited by a harmonic force of amplitude $F_{\text{ext}}$ at the point $(L_x/10, L_y/10)$ and the plate vibratory responses are observed at their centers. In Figure 3, the maximum of the added damping evolves towards a maximum only when $R_L < 1$. The lower the linear resistance in front of 1, the higher the maximum added damping will be for larger values of the relative fluid-solid velocity amplitude. In this configuration, a critical relative velocity is reached, at which the adjusted damping maximum passes through a maximum. For $R_L = 1$, the maximum of the maximum added damping is reached when V tends towards 0. For $R_L > 1$, the maximum added damping is not reached. In this configuration, viscous effects cannot be maximized, whatever the relative fluid-solid velocity.

### 4.3. From a characteristic frequency to a critical relative velocity

In this subsection, the notion of the characteristic frequency provided in Equation (6) is extended to the acoustic nonlinear framework. When $R_L < 1$, the loss factor passes through a maximum at $V_{\text{cr}}$, the critical relative velocity that correspond to the value of V solution of the equation

$$ \frac{\partial \eta_1(\omega, \varepsilon V)}{\partial V} = 0. \qquad (32) $$

The critical relative velocity reads

$$ V_{\text{cr}} = \frac{c_I}{\varepsilon} \left( \frac{\omega}{\omega_c} - 1 \right) \quad \text{with} \quad c_I = \frac{I_1^{(3)}}{I_1^{(4)}} \quad \text{and} \quad \omega \ge \omega_c, \qquad (33) $$
where

$$I_1^{(4)} = \int_S \Psi_1(x, y)^2 |\Psi_1(x, y)| dx dy \qquad (29)$$

carries the Forchheimer's acoustic nonlinear mechanism. As previously mentioned, in addition to being frequency-dependent, the loss factor also depends on the relative fluid-solid velocity in the nonlinear regime. The loss factor provided in Equation (28) is plotted as a function of $\omega$ and $V$ in Figure 3 for a simply-supported MPP of dimensions 490 mm × 570 mm × 1 mm. The perforation ratio, the perforation diameter and the Forchheimer parameter are respectively set to $\phi = 10\%$, $d = 2.2$ mm and $\xi = 0.87 \text{ s m}^{-1}$. In Figure 3, we can see that the maximum of $\eta_1(\omega, V)$ is not unique, but rather corresponds to an infinite

<figure><img src="image_5.png" /><figcaption>Figure 3: Loss factor $\eta_1$ as a function of the angular frequency $\omega$ and the magnitude of the relative fluid-solid velocity $V$. The solid line (—) corresponds to the maximum of $\eta_1$. The simply-supported MPP has dimensions 490 mm × 570 mm × 1 mm, along with a perforation ratio $\phi = 10\%$, a perforation diameter $d = 2.2$ mm, and a Forchheimer parameter $\xi = 0.87 \text{ s m}^{-1}$.</figcaption></figure>

set of points that follow the solid red line. It is important to note here that these observations are only valid with the assumptions made, that is a discretization involving the fundamental harmonic of the first mode only. Under these assumptions, for a fixed $\omega$, there exists a corresponding $V$ that maximizes $\eta_1$, and conversely. The objective is to identify the values of $\omega$ and $V$ that maximize $\eta_1(\omega, V)$, which requires satisfying the condition

$$\mathbf{grad} \eta_1(\omega, V) = 0, \qquad (30)$$

where $\mathbf{grad}(\star) = [\frac{\partial(\star)}{\partial\omega}, \frac{\partial(\star)}{\partial V}]^\top$ is the gradient operator. Equation (30) is solved and after discarding solutions that do not respect the domain of definition of $V$, we obtain

$$\omega_c^{\text{NL}}(V) = \omega_c \left( 1 + \frac{\xi V}{c_I} \right) \quad \text{with} \quad c_I = \frac{I_1^{(3)}}{I_1^{(4)}}, \qquad (31)$$

referred to as the nonlinear characteristic frequency, which extends the concept of the characteristic frequency defined in Equation (6) for the linear regime. Due to the symmetry of $\mathbf{grad} \eta_1(\omega, V)$, Equation (30) admits an infinitely many solutions, and Equation (31) holds true regardless of the value of $V$. This is consistent with the observation made in Figure 3. Equation (31) is rewritten to yield $V$ that satisfies Equation (30) for all $\omega$. The critical velocity reads

$$V_{\text{cr}}(\omega) = \frac{c_I}{\xi} \left( \frac{\omega}{\omega_c} - 1 \right). \qquad (32)$$

In Equations (31) and (32), $\omega_c$ is the characteristic frequency depending on $d$ as defined in Equation (6). Equations (31) and (32) depend on the beam functions arising from the spatial projections of the coupled equations for the first linear plate mode ($I_1^{(3)}$ and $I_1^{(4)}$, encapsulated in $c_I$), the perforation parameters (notably $d$ and $\phi$, included in $\omega_c$), and the Forchheimer parameter $\xi$. As $V$ is defined in Equation (27), $V_{\text{cr}}$ is only defined when $\omega \ge \omega_c$. From the definition of $\omega_c^{\text{NL}}$, the added damping achieves its maximum in the nonlinear regime when $\omega = \omega_c^{\text{NL}}$. For a given MPP, all the parameters in Equation (31) are positive and constant except $V$ which increases with the excitation level. As a result, based on Equation (31), $\omega_c^{\text{NL}}$ can only increase, starting from $\omega_c$, implying that $\omega_c^{\text{NL}} \ge \omega_c$. Three cases can be distinguished according to the value of $\omega$:
<figure><img src="image_5.png" /><figcaption>Figure 3: Loss factor of the first MPP mode obtained by Equation (30) for three MPP: (a) — MPP① corresponding to $R_L > 1$ (—); (b) — MPP② corresponding to $R_L = 1$ (—) and (c) — MPP③ corresponding to $R_L < 1$ (—). The perforation ratio is set to $\phi = 10\%$. The Forchheimer parameter is calculated for each MPP configuration from Equation (9). Corresponding values of $f_c$ obtained from Equation (6) and $R_L$ obtained from Equation (28) as well as plate parameters $(d, L_x, L_y)$ are given in Table 1.</figcaption></figure>

where $\omega_c$ is given by Equation (6). Equation (33) depends on the beam functions in the spatial projections of the coupled equations for the first linear plate mode $I_1^{(3)}$ and $I_1^{(4)}$ and perforation parameters $d$ and $\phi$ in particular. For a specified perforation diameter, the added damping reaches a maximum when $V = V_{cr}$. Therefore, for a given MPP size and knowing $V$, we can adjust $d$ to maximize the added damping.

This behavior has already been observed in the context of sound absorption of an MPP in the nonlinear Forchheimer regime. In fact, the critical relative velocity found in this paper is similar to the critical Mach velocity observed in [7]. The only difference is that in the case presented in this paper, the critical relative velocity given in Equation (33) involves the beam functions in the spatial projections of the coupled equations for the first linear plate mode.

### 4.4. Space distribution of the added damping

Since the added damping depends on the relative fluid-solid velocity, it is therefore a function of space. The present subsection proposes the analysis of the spatial-dependent damping of a simply-supported MPP of dimension 490 mm × 570 mm × 1 mm with $d = 2.2$ mm and $\phi = 10\%$. The single-mode assumption is considered here, and it is assumed that the contributions of harmonics of order greater than one are negligible (c.f Section 3.2). To compute the loss factor and therefore the damping added at each point, the plate is discretized along the $x$ and $y$ axes into $P$ elements of equal size. Each point is defined by its spatial coordinates $(x_g, y_g)$ where $g = 1, 2, \dots, P$. A harmonic excitation of $F_{ext}$ amplitude is applied on the excitation point $(x_e, y_e)$ as given in Equation (12) applied to the considered mode. The Forchheimer parameter is set to $\varepsilon = 0.87 \text{ s m}^{-1}$ from Equation (9). The HBM resolution is performed and outputs the amplitude of the relative fluid-solid velocity at the resonance frequency at each considering point $(x_g, y_g)$. Since the amplitude of the relative fluid-solid velocity is known, the loss factor is obtained from Equation (30).

The overall loss factor at $\omega = \omega_1$ can be estimated by performing the spatial average of the nonlinear loss factor such that

$$ \bar{\eta}_1(\omega_1, \varepsilon V) = \frac{1}{L_x L_y} \int_S \eta_1(\omega_1, \varepsilon V(x, y)) dxdy, \qquad (34) $$

where $\eta_1(\omega_1, \varepsilon V(x, y))$ is given in Equation (30). To this end, a point harmonic excitation of amplitude $F_{ext} \in [10^{-4}; 0.3]$ N is applied and the amplitude of the fluid velocity is calculated for each observation point that corresponds to each point $(x_g, y_g)$ in the plate discretization. Figure 4 plots the normalized overall loss factor of the added damping, given by

$$ L_{\eta_1} = \frac{\bar{\eta}_1(\omega_1, \varepsilon V)}{\eta_1(\omega_1, 0)}, \qquad (35) $$

where $\eta_1(\omega_1, 0)$ represents the loss factor in the linear case obtained for $\varepsilon = 0$, as a function of the spatial mean relative fluid-solid velocity $\langle V \rangle$. The normalized overall loss factor reaches a maximum for a specific value of the mean relative fluid-solid velocity that corresponds to a particular value of the relative fluid-solid amplitude. For small excitations, the overall loss factor of the MPP is about $4.6 \times 10^{-3}$. The
1. $\omega < \omega_c$ — The critical velocity $V_{cr}$ does not exist according to Equation (32). Through Equation (31), $\omega_c^{NL}$ is always greater than $\omega$, regardless of $V$. Consequently, the maximum added damping is reached in neither the linear nor the nonlinear regime, and $\eta_1$ decreases as $V$ increases, as shown below Figure 4(a).

2. $\omega = \omega_c$ — Through Equation (32), the critical velocity $V_{cr}$ vanishes. The maximum added damping occurs when $V$ approaches zero, corresponding to the linear regime. Consequently, $\eta_1(\omega, V)$ reaches its maximum in the linear regime and then decreases with the relative fluid-solid velocity from its linear value, as shown below in Figure 4(b).

3. $\omega > \omega_c$ — The critical velocity is strictly greater than zero according to Equation (32). When the relative velocity $V$ approaches zero, $\omega_c^{NL}$ becomes less than $\omega$. As a result, the maximum viscous dissipation is not achieved in the linear regime but can instead be achieved in the nonlinear regime. Specifically, the added damping for $\omega = \omega_c^{NL}$, does reach its maximum as $V$ increases. Thus, the maximum of $\eta_1$ occurs when $V = V_{cr}$. Consequently, the added damping increases with the relative velocity up to the critical value $V_{cr}$ before decreasing as the relative fluid-solid velocity continues to rise as presented below in Figure 4(c).

To summarize, for a given perforation diameter, the added damping can reach a maximum at $V = V_{cr}$ in the nonlinear regime only if $\omega \ge \omega_c$. Therefore, for specific MPP size and $V$, the diameter $d$ can be adjusted to maximize the added damping. This behavior has already been observed in the context of MPP sound absorption in the nonlinear Forchheimer regime. Notably, the critical relative velocity identified in this study is analogous to the critical Mach velocity observed in [7]<fnref n="5" />. In this paper, the principal difference is that the critical relative velocity given in Equation (32) involves the beam functions in the spatial projections of the coupled equations for the first linear plate mode.

**Remark:** Due to nonlinear damping, the resonance frequency is a function of the forcing magnitude (c.f. Figure 2) and, consequently, of the relative fluid-solid velocity, and is now denoted by $\omega_i^{NL}$ instead of $\omega_i$. For a practical case, the added damping reaches its maximum in

1. the linear regime, when $\omega_i = \omega_c$,

2. the nonlinear regime, when $\omega_i^{NL} = \omega_c^{NL}$.

This is illustrated in Figure 2, where the linear characteristic frequency for the MPP under consideration is $\omega_c = 45 \text{ rad s}^{-1}$, and its first resonance frequency is $\omega_1 = 115 \text{ rad s}^{-1}$. The added damping reaches its maximum around $F_{\text{ext}} = 0.2 \text{ N}$, which corresponds to $V = V_{cr} = 2.44 \text{ m s}^{-1}$ and $\omega_1^{NL} \approx \omega_c^{NL} = 114.1 \text{ rad s}^{-1}$.

In the following, the three cases outlined in the preceding paragraphs are discussed using an example. To this end, the MPP investigated in Section 4.1 is analyzed to illustrate the three cases: $\omega < \omega_c$, $\omega = \omega_c$, and $\omega > \omega_c$. In this analysis, the characteristic frequency in Equation (6) is $\omega_c = 45 \text{ rad s}^{-1}$ and $\omega$ is adjusted to reflect three specific cases: $\omega = 30 \text{ rad s}^{-1}$ for the first case, $\omega = 45 \text{ rad s}^{-1}$ for the second case, and $\omega = 115 \text{ rad s}^{-1}$ for the third case<fnref n="6" />. The perforation ratio and the perforation diameter are respectively set to $\phi = 10\%$ and $d = 2.2 \text{ mm}$. The Forchheimer parameter is obtained from Equation (9). To represent the added damping in the nonlinear regime, the added loss factor is plotted as a function of the amplitude of the relative fluid-solid velocity for the three cases in Figure 4. The studied MPP are excited by a harmonic force with amplitude $F_{\text{ext}}$ at the point $(L_x/10, L_y/10)$, and the plate vibratory responses are observed at their centers. For $\omega < \omega_c$, the maximum of the nonlinear added damping is not reached. In this case, viscous effects cannot be maximized, regardless of the relative fluid-solid velocity (see Figure 4(a)). For $\omega = \omega_c$, the maximum of the nonlinear added damping is reached when $V$ tends towards 0 (see Figure 4(b)). In Figure 4, the added damping evolves towards a maximum only when $\omega > \omega_c$. The lower the ratio $\omega/\omega_c$ is relative to 1, the higher the maximum added damping will be for larger values of the relative fluid-solid velocity amplitude (see Figure 4(c)). In this configuration, a critical relative velocity is reached at which the adjusted damping passes through a maximum.

### 4.3. Space distribution of the added damping

Since the added damping depends on the relative fluid-solid velocity, it is a function of space. The present subsection proposes the analysis of the spatial-dependent damping of a simply-supported MPP of dimensions 490 mm × 570 mm × 1 mm with $d = 2.2 \text{ mm}$ and $\phi = 10\%$. The single-mode assumption is considered here, and it is assumed that the contributions of harmonics of order greater than one are negligible (c.f Section 3.2). To compute the loss factor and therefore the damping added at each point, the plate is discretized along the $x$ and $y$ axes into $P$ elements of equal size. Each point is defined by its spatial coordinates $(x_g, y_g)$ where $g = 1, 2, \dots, P$. A harmonic excitation of $F_{\text{ext}}$ amplitude is applied on

<footnote n="5">⁵For more details, see Equation (22) in [7].</footnote>

<footnote n="6">⁶The latter configuration corresponds to the one illustrated in Figure 2.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4: Normalized overall nonlinear loss factor as a function of the forcing amplitude. The microperforated plate is simplified-supported, and the parameters for the perforations are *d* = 2.2 mm and *φ* = 10%. The Forchheimer parameter is determined using Equation (9), resulting in *ε* = 0.87 s m⁻¹. The excitation amplitudes of points Ⓐ, Ⓑ, Ⓒ and Ⓓ as well as their corresponding coordinates are listed in Table 2.</figcaption></figure>

range in which the normalized overall loss factor remains relatively constant defines the linear regime. In this case, the loss factor obtained with Equation (30) is equal to its value in the linear regime given in Equation (4). It increases with excitation level up to a critical value before decreasing. Finally, Figure 4 shows that there exists an excitation level from which the added damping becomes lower than that obtained in the linear case (for a fixed diameter).

<table><caption>Table 2: Coordinates of interest points of Figure 4 and corresponding external amplitude force F<sub>ext</sub>.</caption><thead><tr><th></th><th>Excitation Ⓐ</th><th>Excitation Ⓑ</th><th>Excitation Ⓒ</th><th>Excitation Ⓓ</th></tr></thead><tbody><tr><td>F<sub>ext</sub> (N)</td><td>0.02</td><td>0.08</td><td>0.15</td><td>0.55</td></tr><tr><td>η<sub>1</sub>(ω<sub>1</sub>, εV)[×10<sup>-3</sup>]</td><td>4.78</td><td>5.13</td><td>5.34</td><td>4.92</td></tr><tr><td>⟨V⟩ (m s<sup>-1</sup>)</td><td>0.09</td><td>0.34</td><td>0.87</td><td>2.21</td></tr></tbody></table>

In Figure 4, four forcing amplitudes of interest are highlighted whose coordinates and corresponding F<sub>ext</sub> are shown in Table 2. For each point of interest, the added loss factor is spatially represented at the first resonance frequency in Figure 5. In Section 4.3, it was shown that the added damping can reach

<figure><img src="image_7.png" /><figcaption>Figure 5: Spatial distribution of the loss factor η₁(ω₁, εV(x, y)) on the plate, obtained from Equation (30) at the first resonance frequency ω₁ for the four points listed in Table 2 and shown in Figure 4.</figcaption></figure>

a maximum at a particular relative fluid-solid velocity. The same results are obtained here. As F<sub>ext</sub> increases, the nonlinearities are activated and the added damping is maximal for the space-dependent critical relative velocity. The higher F<sub>ext</sub>, the greater the relative displacement at the center of the plate. If the critical relative velocity is exceeded, the added damping will decrease at the center of the structure. However, since the velocity is a field, V<sub>cr</sub> is reached at other points of the MPP. The maximum added damping will follow the zone where |ẇ| ≈ V<sub>cr</sub>. This is because damping is governed by the thickness of the boundary layers in the perforations, which are themselves modified when the fluid flow velocity and, therefore, the relative velocity in the case of a vibrating MPP, becomes significant.
<figure><img src="image_6.png" /><figcaption>Figure 4: Loss factor obtained by Equation (28) of the first mode of the simply-supported aluminum MPP of dimension 490 mm × 570 mm × 1 mm with φ = 10 % and d = 2.2 mm. The linear characteristic frequency ω$_{c}$ = 45 rad s − 1 is obtained from Equation (6). Three cases are presented: (a) ω < ω$_{c}$ with ω = 30 rad s − 1 ; (b) ω = ω$_{c}$ and (c) ω > ω$_{c}$ with ω = 115 rad s − 1 . The Forchheimer parameter is calculated from Equation (9).</figcaption></figure>

the excitation point ( x$_{e}$, y$_{e}$ ) as given in Equation (12) applied to the considered mode. The Forchheimer parameter is set to ξ = 0 . 87 s m − 1 from Equation (9). The HBM resolution is performed and outputs the amplitude of the relative fluid-solid velocity at the resonance frequency at each considering point ( x$_{g}$, y$_{g}$ ). Since the amplitude of the relative fluid-solid velocity is known, the loss factor is obtained from Equation (28).

The overall loss factor at ω = ω$_{1}$ can be estimated by performing the spatial average of the nonlinear loss factor such that

$$
\bar{\eta}_1(\omega_1, V) = \frac{1}{L_x L_y} \int_S \eta_1(\omega_1, V(x, y)) \, dx \, dy, \qquad (33)
$$

where η$_{1}$ ( ω$_{1}$, V ( x, y )) is given in Equation (28). To this end, a point harmonic excitation of amplitude F$_{ext}$ ∈ [10 − $^{4}$; 0 . 3] N is applied and the amplitude of the fluid velocity is calculated for each observation point that corresponds to each point ( x$_{g}$, y$_{g}$ ) in the plate discretization. Figure 5 plots the normalized overall loss factor of the added damping, given by

$$
L_{\eta_1} = \frac{\bar{\eta}_1(\omega_1, V)}{\eta_1^L(\omega_1)}, \qquad (34)
$$

where η L 1 represents the loss factor in the linear regime defined in Equation (4). The normalized overall

<figure><img src="image_7.png" /><figcaption>Figure 5: Normalized overall nonlinear loss factor as a function of the forcing amplitude. The microperforated plate is simply-supported, and the parameters for the perforations are d = 2 . 2 mm and φ = 10 %. The Forchheimer parameter is determined using Equation (9), resulting in ξ = 0 . 87 s m − $^{1}$. The excitation amplitudes of points Ⓐ , Ⓑ , Ⓒ and Ⓓ as well as their corresponding coordinates are listed in Table 1.</figcaption></figure>

loss factor is maximum when 〈 V 〉 = V$_{cr}$ . For small excitation, the overall loss factor of the MPP is about 4 . 6 × 10 − $^{3}$. The range in which the normalized overall loss factor remains relatively constant defines the linear regime. In this case, the loss factor obtained with Equation (28) is equal to its value in the linear
# 5. Experimental validation

## 5.1. Experimental set-up

Experimental measurements on samples of free-free microperforated plates excited at the center to reproduce the vibratory behavior of a cantilever beam according to [39], are carried out to validate the analytical model. The experimental set-up is presented in Figure 6. An aluminum microperforated plate

<figure><img src="image_8.png" /><figcaption>Figure 6: Experimental set-up and sample used for the frequency response of microperforated plate: ((a)) and ((b)) represent the experimental set-up respectively in the plan xz and xy; ((c)) provides the microperforated sample investigated experimentally with $x_e$ the abscissa of the excitation point and $x_g$ the abscissa of the observation point. The MPP is made of aluminum and has a dimension of 560 mm × 30 mm × 1 mm. The perforation parameters are set to $d = 2.8$ mm and $\phi = 10\%$.</figcaption></figure>

of dimension 560 mm × 30 mm × 1 mm with $d = 2.8$ mm and $\phi = 10\%$ is considered. The perforation diameter is chosen in order to have $R_L = 0.18 < 1$ and $f_c < f_1$. The plates are chosen to be longer than they are wide, so that the experimentally observed peaks can be easily associated with the proposed theory. In this configuration, mode shapes along $y$ do not participate in the low-frequency range. The previous model is simplified using the Euler-Bernoulli theory. The MPP is excited by a sine-point force along the $z$-axis at the abscissa $x_e = 0$. In this configuration, only even modes of a free-free beam are excited. Its modal behavior is similar to that of a clamped beam of half its length [39]. In fact, the displacement relative to the imposed motion and its spatial derivative are zero at $x_e$

$$w_s(x_e, L_y/2, t) = w(x_e, L_y/2, t) = 0, \quad \text{and} \quad \left. \frac{\partial w_s(x, L_y/2, t)}{\partial x} \right|_{x_e} = \left. \frac{\partial w(x, L_y/2, t)}{\partial x} \right|_{x_e} = 0. \qquad (36)$$

In order to localize the resonance frequencies, the proposed MPP is first excited by white noise to obtain the reference frequencies. Next, the frequency response is measured, which corresponds to the velocity $V_g$ of the solid at point $(x_g, y_g) = (0.26 \text{ m}, L_y/2)$ divided by the base velocity $V_e$ of the solid at the excitation point. To do this, the base acceleration is fixed at the excitation point and must remain constant for the duration of the measurement [40, 41, 42]. During measurements, velocities are measured with a laser Doppler vibrometer (Polytec PSV400 scanning laser vibrometer). Since the excitation is harmonic, the base acceleration $\ddot{w}_e(t)$ is obtained from $2\pi f_e \ddot{w}_e(t)$ where $f_e$ is the forced frequency. The shaker amplitude level is manually adjusted to compensate for and avoid force drops near resonance due to the interaction of the structure with the shaker. The forcing frequency $f_e$ is then adjusted to sweep the frequencies around the first MPP resonance $f_1$ with a frequency step of 0.1 Hz. For each measurement point, 30 measurements are averaged by the laser vibrometer.

## 5.2. Relative displacement

The configuration investigated experimentally corresponds to an MPP subject to an excitation at its base. The displacement measured at the end of the structure is the relative displacement, i.e. the absolute displacement at the end of the beam with respect to the base displacement imposed by the excitation of the point of force. The model gives the total/absolute solid displacement. In order to be able to compare the relative solid motion fields obtained analytically and experimentally, the model presented in Section 2
regime given in Equation (4). It increases with excitation level up to a critical value before decreasing. Figure 5 also indicates that there exists an excitation level above which the added damping is lower than that obtained in the linear case (for a fixed diameter).

<table><caption>Table 1: Coordinates of interest points of Figure 5 and corresponding external amplitude force $F_{ext}$.</caption><thead><tr><th></th><th>Excitation A)</th><th>Excitation B)</th><th>Excitation C)</th><th>Excitation D)</th></tr></thead><tbody><tr><td>F<sub>ext</sub> (N)</td><td>0.02</td><td>0.08</td><td>0.20</td><td>0.55</td></tr><tr><td>η&#773;<sub>1</sub>(ω<sub>1</sub>, V)[×10<sup>-3</sup>]</td><td>4.78</td><td>5.13</td><td>5.34</td><td>4.92</td></tr><tr><td>⟨V⟩ (m s<sup>-1</sup>)</td><td>0.09</td><td>0.34</td><td>0.87</td><td>2.21</td></tr></tbody></table>

In Figure 5, four forcing amplitudes of interest are highlighted whose coordinates and corresponding $F_{ext}$ are shown in Table 1. For each point of interest, the added loss factor is spatially represented at the first resonance frequency in Figure 6. In Section 4.2, it was shown that the added damping can reach

<figure><img src="image_8.png" /><figcaption>Figure 6: Spatial distribution of the loss factor $\eta_1(\omega_1, V(x, y))$ on the plate, obtained from Equation (28) at the first resonance frequency $\omega_1$ for the four excitation points listed in Table 1 and shown in Figure 5: (a) Excitation A) corresponding to $F_{ext} = 0.02$ N; (b) Excitation B) corresponding to $F_{ext} = 0.08$ N; (c) Excitation C) corresponding to $F_{ext} = 0.20$ N; (d) Excitation D) corresponding to $F_{ext} = 0.55$ N.</figcaption></figure>

a maximum at a particular relative fluid-solid velocity. The same results are obtained here. As $F_{ext}$ increases, the nonlinearities are activated and the added damping is maximal for the space-dependent critical relative velocity. The higher $F_{ext}$, the greater the relative displacement at the center of the plate. If the critical relative velocity is exceeded, the added damping will decrease at the center of the structure. However, since the velocity is a field, $V_{cr}$ is reached at other points of the MPP. The maximum added damping will follow the zone where $V = V_{cr}$. This is because damping is governed by the thickness of the boundary layers in the perforations, which are themselves modified when the fluid flow velocity and, therefore, the relative velocity in the case of a vibrating MPP, becomes significant.

# 5. Experimental validation

## 5.1. Experimental set-up

Experimental measurements on samples of free-free microperforated plates excited at the center to reproduce the vibratory behavior of a cantilever beam according to [41], are carried out to validate the analytical model. The experimental set-up is presented in Figure 7. An aluminum microperforated plate of dimensions 560 mm × 30 mm × 1 mm with $d = 2.8$ mm and $\phi = 10\%$ is considered. The perforation diameter is chosen in order to have $f_1 > f_c$. The plates are chosen to be longer than they are wide, so that the experimentally observed peaks can be easily associated with the proposed theory. In this configuration, mode shapes along $y$ do not participate in the low-frequency range. The previous model is simplified using the Euler-Bernoulli theory. The MPP is excited by a sine-point force along the $z$-axis at the abscissa $x_e = 0$. In this configuration, only even modes of a free-free beam are excited. Its modal behavior is similar to that of a clamped beam of half its length [41]. In fact, the displacement relative to the imposed motion and its spatial derivative are zero at $x_e$

$$w_s(x_e, L_y/2, t) = w(x_e, L_y/2, t) = 0, \quad \text{and} \quad \left. \frac{\partial w_s(x, L_y/2, t)}{\partial x} \right|_{x=x_e} = \left. \frac{\partial w(x, L_y/2, t)}{\partial x} \right|_{x=x_e} = 0. \quad (35)$$

In order to identify the resonance frequencies, the proposed MPP is first excited by white noise to obtain the reference frequencies. Next, the frequency response is measured, which corresponds to the velocity $V_g$ of the solid at point $(x_g, y_g) = (0.26 \text{ m}, L_y/2)$ divided by the base velocity $V_e$ of the solid at the excitation point. To do this, the base acceleration is fixed at the excitation point and must remain
is adapted so that the base excitation acts as a forcing point. The total displacement of a structure cross-section can be therefore decomposed as

$$w_s(x, L_y/2, t) = w_e(t) + w_{sr}(x, L_y/2, t) \quad \text{and} \quad w(x, L_y/2, t) = w_e(t) + w_r(x, L_y/2, t), \qquad (37)$$

where $w_e(t)$ is the base displacement. The relative solid displacement and the relative fluid-solid displacement with respect to the base displacement are respectively noted $w_{sr}(x, y, t)$ and $w_r(x, y, t)$. Equation (37) is substituted into Equation (10) which leads to

$$h\rho\ddot{w}_{sr} + h\rho_f\ddot{w}_r + D\nabla^4 w_{sr} = -\ddot{w}_e(h\rho - \phi h\rho_f), \quad (38a)$$

$$\rho_f \ddot{w}_{sr} + \frac{\rho_f \alpha_\infty}{\phi} \ddot{w}_r + \sigma_0 (\dot{w}_r - \phi \dot{w}_e) + \sigma_0 \varepsilon (\dot{w}_r - \phi \dot{w}_e) |(\dot{w}_r - \phi \dot{w}_e)| + \alpha M_f \nabla^2 w_{sr} = -\ddot{w}_e (\rho_f - \rho_f \alpha_\infty). \quad (38b)$$

For a harmonic base excitation, the base acceleration writes

$$\ddot{w}_e(t) = \gamma_e \cos \omega_e t, \qquad (39)$$

where $\omega_e$ is the forcing pulsation, and $\gamma_e$ the amplitude of the base acceleration. The relative displacements $w_{sr}(x, y, t)$ and $w_r(x, y, t)$ are expanded onto the natural mode basis as in Equation (11) and the resulting system of equations is solved following the procedure introduced in Section 3. It should be noted that the nonlinear function in Equation (15) is updated by considering the base displacement, such that

$$f_i^{\text{NL}}(\dot{\mathbf{q}}(t)) = h\sigma_0\varepsilon \int_S \sum_{j=1}^{N} (\dot{q}_j(t)\Psi_j(x, y) + \phi\dot{w}_e(t)) \left| \sum_{p=1}^{N} (\dot{q}_p(t)\Psi_p(x, y) + \phi\dot{w}_e(t)) \right| \Psi_i(x, y)dxdy, \quad (40)$$

where $\dot{\mathbf{q}}(t) = [\dot{q}_1(t), \dot{q}_2(t), \dots, \dot{q}_N(t)]^\top$ stores the temporal contribution of $w_r(x, y, t)$. Experimentally, the resonant response of the structure at the point $x_g$ is measured for several input accelerations and the amplitude at the resonance of the relative solid displacement $w_{sr}(x_g, L_y/2, t)$ is noted $z_r$ and writes from Equation (11)

$$z_r = z_{\max} \frac{\Psi(x_g, L_y/2)}{\Psi_{\max}}, \qquad (41)$$

where $z_{\max}$ is the amplitude of the transversal MPP motion where the deformed shape is maximum. The corresponding value of the mode shape is $\Psi_{\max} = \max|\Psi(x_g, L_y/2)|$. For the first MPP mode of the configuration proposed in Figure 6, the maximum of the deformed shape is reached at the end corner $(L_x, L_y)$ of the structure.

### 5.3. Experimental and analytical results

This subsection presents the experimental results, for non-perforated plate (Figure 7) and for microperforated plate (Figure 6(b)), obtained from the setup presented in Figure 6(a), compared with those derived from the analytical model. The forcing frequency is linked to the base acceleration, which is chosen from the range $\gamma_e \in [10, 90] \text{ mm s}^{-2}$. Firstly, experimental measurement are performed on a non-perforated plate subjected to base excitation. In Figure 7, the normalized velocity $V_g/V_e$ is given as a function of the normalized frequency $f_e/f_1$ with $f_1 = 10.05 \text{ Hz}$. In Figure 6, the base accelerations are set to $\gamma_e = 70 \text{ mm s}^{-2}$ and $\gamma_e = 90 \text{ mm s}^{-2}$. The aim here, is to delimit the range of $\gamma_e$ to be imposed in order to remain in the acoustic nonlinear regime induced only by large fluid displacements. Figure 7 shows that for the $\gamma_e$ values considered, the level amplitude has no influence on the plate transfer function. The vibratory response of the plate without perforation evolves linearly with $\gamma_e$ and therefore with the amplitude of the forcing. Assuming that the structural behavior of the non-perforated plate and the MPP are similar$^6$, it can be stated that for $\gamma_e \le 90 \text{ mm s}^{-2}$, geometric nonlinearities are not activated.

Thus, the experimental measurements for MPP are performed for $\gamma_e$ values lower than $90 \text{ mm s}^{-2}$ in Figure 8. This figure presents the normalized velocity as a function of the normalized frequency, with a comparison between experimental measurement and analytical results for four base accelerations. The forcing frequency is normalized by $f_1 = 9.28 \text{ Hz}$ that corresponds to the resonance frequency with added damping in the linear case. The amplitude variations observed in Figure 8 are therefore only due to the nonlinearities associated with fluid-structure coupling in the perforations. In Figure 8, analytical results are provided by solving Equation (38) using the procedure detailed in Section 3 for $\varepsilon = 1.43 \text{ s m}^{-1}$ which is experimentally measured on a high airflow resistivity meter (c.f. Section 2.3). The value of $\varepsilon$

<footnote>$^6$It was experimentally confirmed that mode shapes were the same between non-perforated plate and microperforated plate.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7: Experimental set-up and sample used for the frequency response of microperforated plate: (a) and (b) represent the experimental set-up respectively in the plan $xz$ and $xy$; (c) provides the microperforated sample investigated experimentally with $x_e$ the abscissa of the excitation point and $x_g$ the abscissa of the observation point. The MPP is made of aluminum and has a dimension of 560 mm × 30 mm × 1 mm. The perforation parameters are set to $d = 2.8$ mm and $\phi = 10\%$.</figcaption></figure>

constant for the duration of the measurement [42–44]. During measurements, velocities are measured with a laser Doppler vibrometer (Polytec PSV400 scanning laser vibrometer). Since the excitation is harmonic, the base acceleration $\ddot{w}_e(t)$ is obtained from $2\pi f_e \dot{w}_e(t)$ where $f_e$ is the forced frequency. The shaker amplitude level is manually adjusted to compensate for and avoid force drops near resonance due to the interaction of the structure with the shaker. The forcing frequency $f_e$ is then adjusted to sweep the frequencies around the first MPP resonance $f_1$ with a frequency step of 0.1 Hz. For each measurement point, 30 measurements are averaged by the laser vibrometer.

## 5.2. Relative displacement

The configuration investigated experimentally corresponds to an MPP subject to an excitation at its base. The displacement measured at the end of the structure is the relative displacement, i.e. the absolute displacement at the end of the beam with respect to the base displacement imposed by the excitation of the point of force. The model gives the total/absolute solid displacement. In order to be able to compare the relative solid motion fields obtained analytically and experimentally, the model presented in Section 2 is adapted so that the base excitation acts as a forcing point. The total displacement of a structure cross-section can therefore be decomposed as

$$w_s(x, L_y/2, t) = w_e(t) + w_{s_r}(x, L_y/2, t) \quad \text{and} \quad w(x, L_y/2, t) = w_e(t) + w_{r}(x, L_y/2, t), \qquad (36)$$

where $w_e(t)$ is the base displacement. The relative solid displacement and the relative fluid-solid displacement with respect to the base displacement are respectively noted $w_{s_r}(x, y, t)$ and $w_r(x, y, t)$. Equation (36) is substituted into Equation (10) which leads to

$$h\rho\ddot{w}_{s_r} + h\rho_f\ddot{w}_r + D\nabla^4 w_{s_r} = -\ddot{w}_e(h\rho - \phi h\rho_f), \quad (37a)$$

$$\rho_f \ddot{w}_{s_r} + \frac{\rho_f \alpha_\infty}{\phi} \ddot{w}_r + \sigma_0 (\dot{w}_r - \phi \dot{w}_e) + \sigma_0 \xi (\dot{w}_r - \phi \dot{w}_e) |(\dot{w}_r - \phi \dot{w}_e)| + \alpha M_f \nabla^2 w_{s_r} = -\ddot{w}_e (\rho_f - \rho_f \alpha_\infty). \quad (37b)$$

For a harmonic base excitation, the base acceleration reads

$$\ddot{w}_e(t) = \gamma_e \cos \omega_e t, \qquad (38)$$

where $\omega_e$ is the forcing pulsation, and $\gamma_e$ the amplitude of the base acceleration. The relative displacements $w_{s_r}(x, y, t)$ and $w_r(x, y, t)$ are expanded onto the natural mode basis as in Equation (11) and the resulting system of equations is solved following the procedure introduced in Section 3. It should be noted that the nonlinear function in Equation (15) is updated by considering the base displacement, such that

$$f_i^{\text{NL}}(\dot{\mathbf{q}}(t)) = h\sigma_0\xi \int_S \sum_{k=1}^{N} (\dot{q}_k(t)\Psi_k(x, y) + \phi\dot{w}_e(t)) \left| \sum_{p=1}^{N} (\dot{q}_p(t)\Psi_p(x, y) + \phi\dot{w}_e(t)) \right| \Psi_i(x, y) \, dx \, dy, \quad (39)$$
<figure><img src="image_9.png" /><figcaption>Figure 7: Non-perforated plate velocity normalized by base velocity around the first plate resonance $f_1$ for two base acceleration (experimental results are represented by points and the analytical result by a solid line): $\gamma_e = 70 \text{ mm s}^{-2}$ (●); $\gamma_e = 90 \text{ mm s}^{-2}$ (▲) and (—). Experimental results are obtained through the set-up presented in Figure 6.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Microperforated plate velocity normalized by base velocity around the first plate resonance $f_1$ for four base acceleration (experimental results shown by points and corresponding analytical results shown by solid or dotted lines): $\gamma_e = 10 \text{ mm s}^{-2}$ (▲); $\gamma_e = 45 \text{ mm s}^{-2}$ (●); $\gamma_e = 60 \text{ mm s}^{-2}$ (●); $\gamma_e = 80 \text{ mm s}^{-2}$ (■). Experimental results obtained through the set-up presented in Figure 6 and solving Equation (38) with $\varepsilon = 1.43 \text{ s m}^{-1}$. The perforation ratio is set to $\phi = 10\%$ and the perforation diameter $d = 2.8 \text{ mm}$ is chosen to have $R_L = 0.18 < 1$.</figcaption></figure>

used corresponds to the average obtained for a series of 8 experimental measurements. The comparison between experimental values and analytical results presented in Figure 8 allows us to validate the proposed model for the range of $\gamma_e$. Their comparison leads to the conclusion that, with regards to the normalized resonance frequency, the stiffness of the system is respectively reduced as the level of the forcing amplitude increases. In the case of the investigated MPP, the perforation diameter is chosen to have $R_L < 1$ with $R_L = 0.18$. As the amplitude of the forcing increases, the normalized amplitude of the plate velocity passes through a minimum that corresponds to a maximum of added damping and $R_{NL} = 1$. Thus, in Figure 6(b), the added damping reaches a maximum around $\gamma_e = 60 \text{ mm s}^{-2}$ before increases with the forcing amplitude.

It can be noticed that the differences between the model and the measurements, noticeable in Figure 8, are probably due to the uncertainty in the measurement of the base acceleration. This is because the displacement amplitude is smaller, and measurement errors are greater.

## 6. Conclusions

This paper proposes a dynamic model to capture the forced vibration behavior of the microperforated plate in the context of the nonlinear acoustic regime, derived from the linear model proposed by the authors in a previous work [5]. The resistive effects due to the high fluid velocity of the airflow in the perforations is captured analytically by using the Forchheimer's correction, a function of the fluid-solid relative velocity in the perforations. Introducing this correction into the governing equation of motion leads to a PDE system with an added damping quadratic term multiplying the relative fluid-solid velocity
where $\dot{\mathbf{q}}(t) = [\dot{q}_1(t), \dot{q}_2(t), \dots, \dot{q}_N(t)]^\top$ stores the temporal contribution of $w_r(x, y, t)$.

Experimentally, the resonant response of the structure at the point $(x_g, L_y/2)$ is measured for several input accelerations and $z_{r_i}$, which corresponding to the amplitude at the $i$th resonance of the relative solid displacement $w_{s_r}(x_g, L_y, t)$, is estimate by [21]:

$$z_{r_i} = z_i^{\max} \frac{\Psi_i(x_g, L_y/2)}{\Psi_i^{\max}}, \qquad (40)$$

where $z_i^{\max}$ is the amplitude of the transverse motion of the structure where the deformed shape for the mode $i$ is maximum and $\Psi_i^{\max} = \max|\Psi_i(x, L_y/2)|$ its corresponding value. For the first beam mode of the configuration proposed in Figure 7, the maximum of the deformed shape is reached at the end $(L_x, L_y/2)$ of the structure.

### 5.3. Experimental and analytical results

This subsection presents the experimental results, for non-perforated plate (Figure 8) and for microp-erforated plate (Figure 7(b)), obtained from the set-up presented in Figure 7(a), compared with those derived from the analytical model. The forcing frequency is linked to the base acceleration, which is chosen from the range $\gamma_e \in [10, 90] \text{ mm s}^{-2}$. Firstly, experimental measurements are performed on a non-perforated plate subjected to base excitation. In Figure 8, the normalized velocity $V_g/V_e$ is given as a function of the normalized frequency $f_e/f_1$ with $f_1 = 10.05 \text{ Hz}$. In Figure 7, the base accelerations are

<figure><img src="image_10.png" /><figcaption>Figure 8: Non-perforated plate velocity normalized by base velocity around the first plate resonance frequency $f_1 = 10.05 \text{ Hz}$, for two base accelerations (experimental results are represented by points and the analytical result by a solid line): $\gamma_e = 70 \text{ mm s}^{-2}$ (●); $\gamma_e = 90 \text{ mm s}^{-2}$ (▲). Experimental results are obtained through the set-up presented in Figure 7.</figcaption></figure>

set to $\gamma_e = 70 \text{ mm s}^{-2}$ and $\gamma_e = 90 \text{ mm s}^{-2}$. The aim here, is to delimit the range of $\gamma_e$ to be imposed in order to remain in the acoustic nonlinear regime induced only by large fluid displacements. Figure 8 shows that for the $\gamma_e$ values considered, the level amplitude has no influence on the plate transfer function. The vibratory response of the plate without perforation evolves linearly with $\gamma_e$ and therefore with the amplitude of the forcing. Assuming that the structural behavior of the non-perforated plate and the MPP are similar$^7$, it can be stated that for $\gamma_e \le 90 \text{ mm s}^{-2}$, geometric nonlinearities are not activated.

Thus, the experimental measurements for MPP are performed for $\gamma_e$ values lower than $90 \text{ mm s}^{-2}$ in Figure 9. This figure presents the normalized velocity as a function of the normalized frequency, with a comparison between experimental measurement and analytical results for four base accelerations. The forcing frequency is normalized by $f_1 = 9.28 \text{ Hz}$ that corresponds to the resonance frequency with added damping in the linear case. The amplitude variations observed in Figure 9 are therefore only due to the nonlinearities associated with fluid-structure coupling in the perforations. In Figure 9, analytical results are provided by solving Equation (37) using the procedure detailed in Section 3 for $\xi = 1.43 \text{ s m}^{-1}$ which is experimentally measured on a high airflow resistivity meter (c.f. Section 2.3). The value of $\xi$ used corresponds to the average obtained for a series of 8 experimental measurements. The comparison between experimental values and analytical results presented in Figure 9 allows us to validate the proposed model for the range of $\gamma_e$. Their comparison leads to the conclusion that, with regards to the normalized resonance frequency, the stiffness of the system is respectively reduced as the level of the forcing amplitude

<footnote>$^7$It was experimentally confirmed that mode shapes of the non-perforated and microp perforated plate are the (almost) same.</footnote>
in the equation of motion associated with fluid-structure coupling. The resulting system is numerically solved in the forced regime using a harmonic method.

For an MPP the maximum added damping is for when a resonance frequency coincides with the characteristic frequency. The value of this maximum damping varies as a function of the relative fluid-plate velocity. It has been shown that the maximum viscothermal effect is reached when the normalized nonlinear resistivity is equal to 1. Thus, if a resonance frequency coincides with the characteristic frequency and the normalized nonlinear resistivity is equal to 1, a characteristic relative fluid-solid velocity is defined for which the added damping reaches a maximum. Since the normalized nonlinear resistance is a linear increasing function of relative speed, there are three categories of MPP : (i) MPP with a normalized linear resistance of less than 1, in which case the maximum damping is reached in the nonlinear regime, (ii) MPP with a normalized linear resistance equal to 1, in which case, the maximum of the damping maximum is reached in the linear regime and decreased with the relative fluid-solid velocity, and (iii) MPP with a normalized linear resistance greater than one, in which case the viscous dissipative effects are not maximized in the nonlinear regime.

As the relative fluid-solid velocity is a spatial function, the maximum nonlinear added damping follows the zone where the relative velocity is equal to the characteristic relative velocity. It is found a value of the forcing amplitude for which the overall averaging nonlinear damping is maximum.

The proposed model is validated by measurements on a microperforated plate.

As a final comment, it should be mentioned that other phenomena must be considered making an accurate model for higher levels of excitation. First, in the case of large deflections, two other types of non-linearity must be added to the model: (i) geometric non-linearity, to reflect the fact that stresses are no longer proportional to deformations, and (ii) quadratic drag damping, to account for friction with surrounding air on both sides of the structure's surface.

## CRediT authorship contribution statement

**Lucie Gallerand:** Conceptualization, Methodology, Software, Validation, Writing - original draft, Writing - review & editing. **Mathias Legrand:** Methodology, Software, Supervision, Writing - original draft, review & editing. **Raymond Panneton:** Writing - original draft, review & editing. **Philippe Leclaire:** Conceptualization, Methodology, Supervision, Writing - original draft, review & editing. **Thomas Dupont:** Conceptualization, Methodology, Supervision, Resources, Funding acquisition, Writing - original draft, review & editing.

## Acknowledgment

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the project of the reference number RGPIN-2019-06573 and by Fonds de recherche du Québec - Nature et technologies (FRQNT) program financing.

## References

[1] M.Q. Wu. Micro-perforated panels for duct silencing. *Noise Control Engineering Journal*, 45(2):69-77, 1997. [DOI].

[2] L. Wang and T. Wang. Investigation of the effect of perforated sheath on thermal-flow characteristics over a gas turbine reverse-flow combustor—Part 2: Computational analysis. *Journal of Thermal Science and Engineering Applications*, 12(4), 2019. [DOI], [OA].

[3] D.-Y. Maa. Potential of microperforated panel absorber. *Journal of the Acoustical Society of America*, 104(5):2861-2866, 1997. [DOI], [OA].

[4] N. Atalla and F. Sgard. Modeling of perforated plates and screens using rigid frame porous models. *Journal of Sound and Vibration*, 303:195-208, 2007. [DOI], [OA].

[5] L. Gallerand, M. Legrand, T. Dupont, and P. Leclaire. Vibration and damping analysis of a thin finite-size microperforated plate. *Journal of Sound and Vibration*, 541:117295-1-20, 2022. [DOI], [OA].

[6] L. Gallerand, M. Legrand, T. Dupont, R. Panneton, and P. Leclaire. Damping performance of finite microperforated plates using multi-size and spatial distributions of perforations. *Applied Acoustics*, 2024. [OA].

[7] R. Tayong, T. Dupont, and P. Leclaire. On the variations of acoustic absorption peak with particle velocity in micro-perforated panels at high level of excitation. *The Journal of the Acoustical Society of America*, 127(5):2875-2882, 2010. [DOI], [OA].

[8] M. A. Temiz, J. Tournadre, I. L. Arteaga, and A. Hirschberg. Non-linear acoustic transfer impedance of micro-perforated plates with circular orifices. *Journal of Sound and Vibration*, 366:418-428, 2016. [DOI], [OA].

[9] D.-Y. Maa. Microperforated panel at high sound intensity. In *INTER-NOISE and NOISE-CON Congress and Conference*, volume 1994, pages 1511-1514. Institute of Noise Control Engineering, 1994. [OA].

[10] A. Cummings. Acoustic nonlinearities and power losses at orifices. *AIAA Journal*, 22(6):786-792, 1984. [DOI].
<figure><img src="image_11.png" /><figcaption>Figure 9: Microperforated plate velocity normalized by base velocity around the first plate resonance frequency $f_1 = 9.28$ Hz, for four base accelerations (experimental results shown by points and corresponding analytical results shown by solid or dotted lines): $\gamma_e = 10 \text{ mm s}^{-2}$ (▲); $\gamma_e = 45 \text{ mm s}^{-2}$ (●); $\gamma_e = 60 \text{ mm s}^{-2}$ (■); $\gamma_e = 80 \text{ mm s}^{-2}$ (■). Experimental results obtained through the set-up presented in Figure 7 and solving Equation (37) with $\xi = 1.43 \text{ s m}^{-1}$. The perforation ratio is $\phi = 10\%$ and the perforation diameter $d = 2.8$ mm is chosen to have $f_1 > f_c$.</figcaption></figure>

increases. In the case of the investigated MPP, the perforation diameter is chosen to have $f_c < f_1$. As the amplitude of the forcing increases, the normalized amplitude of the plate velocity passes through a minimum that corresponds to a maximum of added damping. Thus, in Figure 9, the added damping reaches a maximum around $\gamma_e = 60 \text{ mm s}^{-2}$.

It can be noticed that the differences between the model and the measurements, noticeable in Figure 9, are probably due to the uncertainty in the measurement of the base acceleration. This is because the displacement amplitude is smaller, and measurement errors are greater.

## 6. Conclusions

This paper presents a dynamic model designed to capture the forced response of a microperforated plate operating within a nonlinear acoustic regime, building upon the linear model previously developed by the authors [5]. The resistive effects due to the high fluid velocity of the airflow in the perforations is introduced analytically by the Forchheimer correction, which is a function of the fluid-solid relative velocity in the perforations. The resulting equation of motion bears an added damping quadratic term multiplying the relative fluid-solid velocity associated with fluid-structure coupling. The system is numerically solved in the forced regime using the harmonic balance method.

It is found that the maximum added damping is reached when a resonance frequency of the system coincides with the characteristic frequency, derived from Biot's frequency for porous materials. This maximum damping varies as a function of the relative fluid-solid velocity. In the present work, the emphasis is placed on the resonance frequency. It is demonstrated that the maximum added damping, driven by the viscothermal dissipation mechanism, can be achieved in the nonlinear regime when the considering driving frequency exceeds the characteristic frequency. This maximum occurs at a critical relative fluid-solid velocity, determined by the perforation characteristics, driving frequency and plate parameters. In a practical case, MPP can be categorized into three groups based on their characteristic frequency: i) characteristic frequency lower than the first resonance frequency; the maximum damping is achieved in the nonlinear regime, ii) characteristic frequency equal to the first resonance frequency; the maximum damping is attained in the linear regime and decreases as the relative fluid-solid velocity increases, and iii) characteristic frequency higher than the first resonance frequency; the added damping does not reach its maximum. As the relative fluid-solid velocity is a spatial function, the maximum nonlinear added damping follows the zone where the relative velocity is equal to the critical velocity. The proposed model is validated by measurements.

As a final comment, it should be mentioned that, in addition to the Forchheimer correction, other mechanisms must be considered to make an accurate model for higher levels of excitation. First, in the case of large deflections, two other types of non-linearity must be added to the model: i) geometric non-linearity, to reflect the fact that stresses are no longer proportional to deformations, and ii) quadratic drag damping, to account for friction with surrounding air on both sides of the structure's surface.
[11] A. Cummings. Transient and multiple frequency sound transmission through perforated plates at high amplitude. *Journal of the Acoustical Society of America*, 79(4):942–951, 1986. [DOI].

[12] T.H. Melling. The acoustic impedance of perforates at medium and high sound pressure levels. *Journal of Sound and Vibration*, 29(1):1–65, 1973. [DOI].

[13] U. Ingard. Absorption Characteristics of Nonlinear Acoustic Resonators. *The Journal of the Acoustical Society of America*, 44(4):1155–1156, 1968. [DOI].

[14] P. Forchheimer. Wasserbewegung durch boden. *Zeitschrift des Vereines Deutscher Ingenieure*, 45(50):1781–1788, 1901.

[15] M. Firdaouss, J.-L. Guermond, and P. Le Quéré. Nonlinear corrections to Darcy’s law at low Reynolds numbers. *Journal of Fluid Mechanics*, 343:331–350, 1997. [DOI].

[16] P. M. Adler, A. E. Malevich, and V. V. Mityushev. Nonlinear correction to Darcy’s law for channels with wavy walls. *Acta Mechanica*, 224(8):1823–1848, 2013. [DOI].

[17] R. Tayong, T. Dupont, and P. Leclaire. Experimental investigation of holes interaction effect on the sound absorption coefficient of micro-perforated panels under high and medium sound levels. *Applied Acoustics*, 72(10):777–784, 2011. [DOI], [OA].

[18] Z. Laly, N. Atalla, and S-A. Meslioui. Acoustical modeling of micro-perforated panel at high sound pressure levels using equivalent fluid approach. *Journal of Sound and Vibration*, 427:134–158, 2018. [DOI].

[19] A. Veisiara, H. Mohammad-Sedighi, and A. Reza. Computational analysis of the nonlinear vibrational behavior of perforated plates with initial imperfection using NURBS-based isogeometric approach. *Journal of Computational Design and Engineering*, 8(5):1307–1331, 2021. [DOI].

[20] M. Amabili. Nonlinear vibrations of rectangular plates with different boundary conditions: theory and experiments. *Computers & Structures*, 82(31-32):2587–2605, 2004. [DOI].

[21] M. Colin, O. Thomas, S. Grondel, and E. Cattan. Very large amplitude vibrations of flexible structures: Experimental identification and validation of a quadratic drag damping model. *Journal of Fluids and Structures*, 97:103056, 2020. [DOI], [OA].

[22] M. Ducceschi, C. Touzé, S. Bilbao, and C. J. Webb. Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations. *Acta Mechanica*, 225(1):213–232, 2014. [DOI], [OA].

[23] S. Bilbao, O. Thomas, C. Touzé, and M. Ducceschi. Conservative numerical methods for the full von Kármán plate equations. *Numerical Methods for Partial Differential Equations*, 31(6):1948–1970, 2015. [DOI], [OA].

[24] V. Denis, A. Pelat, C. Touzé, and F. Gautier. Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity. *International Journal of Non-Linear Mechanics*, 94:134–145, 2017. [DOI] [OA].

[25] H. Li, C. Touzé, F. Gautier, and A. Pelat. Linear and nonlinear dynamics of a plate with acoustic black hole, geometric and contact nonlinearity for vibration mitigation. *Journal of Sound and Vibration*, 508:116206, 2021. [DOI], [OA].

[26] P. Leclaire, K.V. Horoshenkov, and A. Cummings. Transverse vibration of a thin rectangular porous plate saturated by a fluid. *Journal of Sound and Vibration*, 247(1):1–18, 2001. [DOI], [OA].

[27] M.A. Biot and D.G. Willis. The elastic coefficients of the theory of consolidation. *Journal of Applied Mechanics*, pages 594–604, 1957. [DOI].

[28] A.R. Boccaccini and Z. Fan. A new approach for the Young’s modulus-porosity correlation of ceramic materials. *Ceramics International*, 23, 1997. [DOI], [OA].

[29] J.C. Wang. Young’s modulus of porous materials: Part 1 theoretical derivation of modulus-porosity correlation. *Journal of Materials Science*, 19:801–808, 1984. [DOI].

[30] M.A. Biot. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. *Journal of the Acoustical Society of America*, 28(2):168–178, 1956. [DOI], [OA].

[31] Y. Aurégan and M. Pachebat. Measurement of the nonlinear behavior of acoustical rigid porous materials. *Physics of Fluids*, 11(6):1342–1345, 1999. [DOI], [OA].

[32] O. Umnova, K. Attenborough, E. Standley, and A. Cummings. Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment. *The Journal of the Acoustical Society of America*, 114(3), 2003. [DOI].

[33] D. C. Brooke, O. Umnova, P. Leclaire, and T. Dupont. Acoustic metamaterial for low frequency sound absorption in linear and nonlinear regimes. *Journal of Sound and Vibration*, 485, 2020. [DOI].

[34] L. Liu, J.P. Thomas, E.H. Dowell, P. Attar, and K.C. Hall. A comparison of classical and high dimensional harmonic balance approaches for a Duffing oscillator. *Journal of Computational Physics*, 215(1):298–320, 2006. [DOI].

[35] Z.K. Peng, G. Meng, Z.Q. Lang, W.M. Zhang, and F.L. Chu. Study of the effects of cubic nonlinear damping on vibration isolations using Harmonic Balance Method. *International Journal of Non-Linear Mechanics*, 47(10):1073–1080, 2012. [DOI].

[36] T. Zhu, G. Zhang, and C. Zang. Frequency-domain nonlinear model updating based on analytical sensitivity and the Multi-Harmonic balance method. *Mechanical Systems and Signal Processing*, 163:108169, 2022. [DOI].

[37] M. Berthillier, C. Dupont, R. Mondal, and J. J. Barrau. Blades forced response analysis with friction dampers. *Journal of Vibration and Acoustics*, 120(2):468–474, 1998. [DOI], [OA].

[38] E. Pescheck and C. Pierre. An analysis of a friction damped system using two component mode methods. In *International Design Engineering Technical Conferences and Computers and Information in Engineering Conference*, volume 1C: 16th Biennial Conference on Mechanical Vibration and Noise, pages 1–11, 1997. [DOI], [OA].
## CRediT authorship contribution statement

**Lucie Gallerand:** Conceptualization, Methodology, Software, Validation, Writing - original draft, Writing - review & editing. **Mathias Legrand:** Methodology, Software, Supervision, Writing - original draft, review & editing. **Raymond Panneton:** Writing - original draft, review & editing. **Philippe Leclaire:** Conceptualization, Methodology, Supervision, Writing - original draft, review & editing. **Thomas Dupont:** Conceptualization, Methodology, Supervision, Resources, Funding acquisition, Writing - original draft, review & editing.

## Supplementary material

Python scripts used to perform the modal analysis presented in this article are available in commit 004d8318 of the [Git repository] Acoustic nonlinear

## Acknowledgment

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the project of the reference number RGPIN-2019-06573 and by Fonds de recherche du Québec - Nature et technologies (FRQNT) program financing.

## References

[1] M. Wu, Micro-perforated panels for duct silencing, Noise Control Engineering Journal 45 (2) (1997) 69-77,
DOI:10.3397/1.2828428.

[2] L. Wang, T. Wang, Investigation of the effect of perforated sheath on thermal-flow characteristics over a gas turbine reverse-flow combustor—Part 2: Computational analysis, Journal of Thermal Science and Engineering Applications 12 (4) (2019) 041004-1-11, DOI:10.1115/1.4045391.

[3] D.-Y. Maa, Potential of microperforated panel absorber, Journal of the Acoustical Society of America 104 (5) (1997) 2861-2866, DOI:10.1121/1.423870.

[4] N. Atalla, F. Sgard, Modeling of perforated plates and screens using rigid frame porous models, Journal of Sound and Vibration 303 (2007) 195-208, DOI:10.1016/j.jsv.2007.01.012.

[5] L. Gallerand, M. Legrand, T. Dupont, P. Leclaire, Vibration and damping analysis of a thin finite-size microperforated plate, Journal of Sound and Vibration 541 (2022) 117295-1-20, DOI:10.1016/j.jsv.2022.117295.

[6] L. Gallerand, M. Legrand, T. Dupont, R. Panneton, P. Leclaire, Damping performance of finite microperforated plates using multi-size and spatial distributions of perforations, Applied Acoustics 222 (2024) 110041, DOI:10.1016/j.apacoust.2024.110041.

[7] R. Tayong, T. Dupont, P. Leclaire, On the variations of acoustic absorption peak with particle velocity in micro-perforated panels at high level of excitation, The Journal of the Acoustical Society of America 127 (5) (2010) 2875-2882, DOI:10.1121/1.3372714.

[8] M. A. Temiz, J. Tournadre, I. L. Arteaga, A. Hirschberg, Non-linear acoustic transfer impedance of micro-perforated plates with circular orifices, Journal of Sound and Vibration 366 (2016) 418-428, DOI:10.1016/j.jsv.2015.12.022.

[9] D.-Y. Maa, Microperforated panel at high sound intensity, in: INTER-NOISE and NOISE-CON Congress and Conference, Vol. 1994, Institute of Noise Control Engineering, 1994, pp. 1511-1514, OAI:hal-04551854.

[10] A. Cummings, Acoustic nonlinearities and power losses at orifices, AIAA Journal 22 (6) (1984) 786-792,
DOI:10.2514/3.8680.

[11] A. Cummings, Transient and multiple frequency sound transmission through perforated plates at high amplitude, Journal of the Acoustical Society of America 79 (4) (1986) 942-951, DOI:10.1121/1.393691.

[12] T. Melling, The acoustic impedance of perforates at medium and high sound pressure levels, Journal of Sound and Vibration 29 (1) (1973) 1-65, DOI:10.1016/S0022-460X(73)80125-7.

[13] U. Ingard, Absorption Characteristics of Nonlinear Acoustic Resonators, The Journal of the Acoustical Society of America 44 (4) (1968) 1155-1156, DOI:10.1121/1.1911215.

[14] P. Forchheimer, Wasserbewegung durch boden, Zeitschrift des Vereines Deutscher Ingenieure 45 (50) (1901) 1781-1788.

[15] M. Firdaouss, J.-L. Guermond, P. Le Quéré, Nonlinear corrections to Darcy's law at low Reynolds numbers, Journal of Fluid Mechanics 343 (1997) 331-350, DOI:10.1017/S0022112097005843.

[16] P. M. Adler, A. E. Malevich, V. V. Mityushev, Nonlinear correction to Darcy's law for channels with wavy walls, Acta Mechanica 224 (8) (2013) 1823-1848, DOI:10.1007/s00707-013-0840-3.

[17] R. Tayong, T. Dupont, P. Leclaire, Experimental investigation of holes interaction effect on the sound absorption coefficient of micro-perforated panels under high and medium sound levels, Applied Acoustics 72 (10) (2011) 777-784, DOI:10.1016/j.apacoust.2011.04.011.

[18] Z. Laly, N. Atalla, S.-A. Meslioui, Acoustical modeling of micro-perforated panel at high sound pressure levels using equivalent fluid approach, Journal of Sound and Vibration 427 (2018) 134-158, DOI:10.1016/j.jsv.2017.09.011.
[39] J.-L. Wojtowicki, L. Jaouen, and R. Panneton. New approach for the measurement of damping properties of materials using the oberst beam. *Review of scientific instruments*, 75(8):2569–2574, 2004. [DOI], [OA].

[40] G. Zhang, X. Wang, and Z. Yang. Study on excitation force characteristics in a coupled shaker-structure system considering structure modes coupling. *Chinese Journal of Aeronautics*, 35(7):227–245, 2022. [DOI].

[41] B. R. Pacini, R. J. Kuether, and D. R. Roettgen. Shaker-structure interaction modeling and analysis for nonlinear force appropriation testing. *Mechanical Systems and Signal Processing*, 162:108000, 2022. [DOI], [OA].

[42] P. S. Varoto and L. P. R. De Oliveira. On the force drop off phenomenon in shaker testing in experimental modal analysis. *Shock and Vibration*, 9(4-5):165–175, 2002. [DOI].
[19] A. Veisiara, H. Mohammad-Sedighi, A. Reza, Computational analysis of the nonlinear vibrational behavior of perforated plates with initial imperfection using NURBS-based isogeometric approach, Journal of Computational Design and Engineering 8 (5) (2021) 1307-1331, doi:10.1093/jcde/qwab043.

[20] M. Amabili, Nonlinear vibrations of rectangular plates with different boundary conditions: theory and experiments, Computers & Structures 82 (31-32) (2004) 2587-2605, doi:10.1016/j.compstruc.2004.03.077.

[21] M. Colin, O. Thomas, S. Grondel, E. Cattan, Very large amplitude vibrations of flexible structures: Experimental identification and validation of a quadratic drag damping model, Journal of Fluids and Structures 97 (2020) 103056, doi:10.1016/j.jfluidstructs.2020.103056.

[22] M. Ducceschi, C. Touzé, S. Bilbao, C. J. Webb, Nonlinear dynamics of rectangular plates: investigation of modal interaction in free and forced vibrations, Acta Mechanica 225 (1) (2014) 213-232, doi:10.1007/s00707-013-0931-1.

[23] S. Bilbao, O. Thomas, C. Touzé, M. Ducceschi, Conservative numerical methods for the full von Kármán plate equations, Numerical Methods for Partial Differential Equations 31 (6) (2015) 1948-1970, doi:10.1002/num.21974.

[24] V. Denis, A. Pelat, C. Touzé, F. Gautier, Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity, International Journal of Non-Linear Mechanics 94 (2017) 134-145, doi:10.1016/j.ijnonlinmec.2016.11.012.

[25] H. Li, C. Touzé, F. Gautier, A. Pelat, Linear and nonlinear dynamics of a plate with acoustic black hole, geometric and contact nonlinearity for vibration mitigation, Journal of Sound and Vibration 508 (2021) 116206, doi:10.1016/j.jsv.2021.116206.

[26] P. Leclaire, K. Horoshenkov, A. Cummings, Transverse vibration of a thin rectangular porous plate saturated by a fluid, Journal of Sound and Vibration 247 (1) (2001) 1-18, doi:10.1006/jsvi.2001.3656.

[27] M. Biot, D. Willis, The elastic coefficients of the theory of consolidation, Journal of Applied Mechanics (1957) 594-604, doi:10.1115/1.4011606.

[28] A. Boccaccini, Z. Fan, A new approach for the Young's modulus-porosity correlation of ceramic materials, Ceramics International 23, doi:10.1016/S0272-8842(96)00033-8 (1997).

[29] J. Wang, Young's modulus of porous materials: Part 1 theoretical derivation of modulus-porosity correlation, Journal of Materials Science 19 (1984) 801-808, doi:10.1007/BF00540451.

[30] M. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, Journal of the Acoustical Society of America 28 (2) (1956) 168-178, doi:10.1121/1.1908239.

[31] Y. Aurégan, M. Pachebat, Measurement of the nonlinear behavior of acoustical rigid porous materials, Physics of Fluids 11 (6) (1999) 1342-1345, doi:10.1063/1.869999.

[32] O. Umnova, K. Attenborough, E. Standley, A. Cummings, Behavior of rigid-porous layers at high levels of continuous acoustic excitation: Theory and experiment, The Journal of the Acoustical Society of America 114 (3), doi:10.1121/1.1603236 (2003).

[33] D. C. Brooke, O. Umnova, P. Leclaire, T. Dupont, Acoustic metamaterial for low frequency sound absorption in linear and nonlinear regimes, Journal of Sound and Vibration 485, doi:10.1016/j.jsv.2020.115585 (2020).

[34] L. Liu, J. Thomas, E. Dowell, P. Attar, K. Hall, A comparison of classical and high dimensional harmonic balance approaches for a Duffing oscillator, Journal of Computational Physics 215 (1) (2006) 298-320, doi:10.1016/j.jcp.2005.10.026.

[35] Z. Peng, G. Meng, Z. Lang, W. Zhang, F. Chu, Study of the effects of cubic nonlinear damping on vibration isolations using Harmonic Balance Method, International Journal of Non-Linear Mechanics 47 (10) (2012) 1073-1080, doi:10.1016/j.ijnonlinmec.2011.09.013.

[36] T. Zhu, G. Zhang, C. Zang, Frequency-domain nonlinear model updating based on analytical sensitivity and the Multi-Harmonic balance method, Mechanical Systems and Signal Processing 163 (2022) 108169, doi:10.1016/j.ymssp.2021.108169.

[37] M. Berthillier, C. Dupont, R. Mondal, J. J. Barrau, Blades forced response analysis with friction dampers, Journal of Vibration and Acoustics 120 (2) (1998) 468-474, doi:10.1115/1.2893853.

[38] E. Pescheck, C. Pierre, An analysis of a friction damped system using two component mode methods, in: International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. 1C: 16th Biennial Conference on Mechanical Vibration and Noise, 1997, pp. 1-11, doi:10.1115/DETC97/VIB-3844.

[39] O. Thomas, J. Ducarne, J.-F. Deü, Performance of piezoelectric shunts for vibration reduction, Smart Materials and Structures 21 (1) (2011) 015008, https://dx.doi.org/10.1088/0964-1726/21/1/015008.

[40] M. Vakilinejad, A. Grolet, O. Thomas, A comparison of robustness and performance of linear and nonlinear Lanchester dampers, Nonlinear Dynamics 100 (1) (2020) 269-287, doi:10.1007/s11071-020-05512-x.

[41] J.-L. Wojtowicki, L. Jaouen, R. Panneton, New approach for the measurement of damping properties of materials using the oberst beam, Review of scientific instruments 75 (8) (2004) 2569-2574, doi:10.1063/1.1777382.

[42] G. Zhang, X. Wang, Z. Yang, Study on excitation force characteristics in a coupled shaker-structure system considering structure modes coupling, Chinese Journal of Aeronautics 35 (7) (2022) 227-245, doi:10.1016/j.cja.2021.11.024.

[43] B. Pacini, R. Kuether, D. Roettgen, Shaker-structure interaction modeling and analysis for nonlinear force appropriation testing, Mechanical Systems and Signal Processing 162 (2022) 108000, doi:10.1016/j.ymssp.2021.108000.

[44] P. S. Varoto, L. P. R. De Oliveira, On the force drop off phenomenon in shaker testing in experimental modal analysis, Shock and Vibration 9 (4-5) (2002) 165-175, doi:10.1155/2002/675674.