# Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions

Antoine Rallu, Claude Boutin

► To cite this version:

Antoine Rallu, Claude Boutin. Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions. Mathematics, 2024, 12 (12), pp.1844. 10.3390/math12121844. hal-04625202

HAL Id: hal-04625202

https://hal.science/hal-04625202

Submitted on 26 Jun 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions

Antoine Rallu, Claude Boutin

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

## To cite this version:

Antoine Rallu, Claude Boutin. Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions. Mathematics, 2024, 12 (12), pp.1844. 10.3390/math12121844. hal-04625202

HAL Id: hal-04625202

https://hal.science/hal-04625202v1

Submitted on 26 Jun 2024

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
Article

# Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions

Antoine Rallu* and Claude Boutin

Ecole Nationale des Travaux Publics de l'Etat, Université de Lyon, LTDS UMR CNRS 5513, Rue Maurice Audin, 69518 Vaulx-en-Velin, France; claude.boutin@entpe.fr

* Correspondence: antoine.rallu@entpe.fr

**Abstract:** This paper deals with the long-wavelength behaviour of a Euler beam periodically supported by co-located rotation and compression springs. An asymptotic homogenization method is applied to derive the several macroscopic models according to the stiffness contrasts between the elastic supports and the beam. Effective models of differential order two or four are obtained, which can be merged into a single unified model whose dispersion relations at long and medium wavelengths fit those derived by Floquet-Bloch. Moreover, the essential role of rotation supports is clearly evidenced. A mixed "discrete/continuous" approach to the boundary conditions is proposed, which allows the boundary conditions actually applied at the local scale to be expressed in terms of Robin-type boundary conditions on macroscopic variables. This approach can be applied to both dominant-order and higher-order models. The modal analysis performed with these boundary conditions and the homogenised models gives results in good agreement with a full finite element calculation, with great economy of numerical resources.

**Keywords:** homogenisation; dynamics; boundary conditions; asymptotic method; multi-supported beams

MSC: 35B27;74H10;74Q10;74Q15

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

check for
updates

**Citation:** Rallu, A.; Boutin, C. Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions. Mathematics **2024**, *12*, 1844. https://doi.org/10.3390/math12121844

Academic Editors: Julius Kaplunov, Igor Andrianov and Nikolaos L. Tsitsas

Received: 16 April 2024

Revised: 4 June 2024

Accepted: 7 June 2024

Published: 13 June 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Multi-supported beams are widely used in a variety of engineering fields, including long multi-supported bridges, multi-span steel structures, rails connected to sleepers, pipelines supported by regularly distributed foundations, and pipe networks attached to structural elements. Understanding their vibratory behaviour is of interest for safer and more efficient design, in particular for better control of their response to moving loads or to stresses induced by external actions such as ambient vibrations, earthquakes, etc. Another field of application is the monitoring and control of the in-service condition of these structures using structural health monitoring (SHM) vibration methods.

Due to their wide range of applications, many studies have focused on dynamics of beams that rest on continuous foundations or on discrete flexible supports. For good reviews in this field one can refer to [1,2].

As for continuous foundations, the precursor works of [3–5] deal with Euler or Timoshenko beams, respectively, lying on an elastic or visco-elastic foundation of a Winkler compressional type. Let us mention the instability analysis developed by [6] and the work of [7,8] dealing with non-linear foundations and stratified foundations, respectively. The continuous approach has also been applied to discrete foundations by redistributing discrete forces in continuous form, e.g., [9].

Studies that explicitly take into account the discrete nature of supports were initiated by Mead for periodic [10] or non-periodic [11] distributions with the aim of determining their vibro-acoustic behaviour. More recently, ref. [2] investigates the steady-state dynamic
Article

# Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions

Antoine Rallu* and Claude Boutin

Ecole Nationale des Travaux Publics de l'Etat, Université de Lyon, LTDS UMR CNRS 5513, Rue Maurice Audin, 69518 Vaulx-en-Velin, France; claude.boutin@entpe.fr

* Correspondence: antoine.rallu@entpe.fr

**Abstract:** This paper deals with the long-wavelength behaviour of a Euler beam periodically supported by co-located rotation and compression springs. An asymptotic homogenization method is applied to derive the several macroscopic models according to the stiffness contrasts between the elastic supports and the beam. Effective models of differential order two or four are obtained, which can be merged into a single unified model whose dispersion relations at long and medium wavelengths fit those derived by Floquet-Bloch. Moreover, the essential role of rotation supports is clearly evidenced. A mixed "discrete/continuous" approach to the boundary conditions is proposed, which allows the boundary conditions actually applied at the local scale to be expressed in terms of Robin-type boundary conditions on macroscopic variables. This approach can be applied to both dominant-order and higher-order models. The modal analysis performed with these boundary conditions and the homogenised models gives results in good agreement with a full finite element calculation, with great economy of numerical resources.

**Keywords:** homogenisation; dynamics; boundary conditions; asymptotic method; multi-supported beams

MSC: 35B27;74H10;74Q10;74Q15

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

check for
updates

**Citation:** Rallu, A.; Boutin, C. Modal Analysis of a Multi-Supported Beam: Macroscopic Models and Boundary Conditions. Mathematics **2024**, *12*, 1844. https://doi.org/10.3390/math12121844

Academic Editors: Julius Kaplunov, Igor Andrianov and Nikolaos L. Tsitsas

Received: 16 April 2024

Revised: 4 June 2024

Accepted: 7 June 2024

Published: 13 June 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Multi-supported beams are widely used in a variety of engineering fields, including long multi-supported bridges, multi-span steel structures, rails connected to sleepers, pipelines supported by regularly distributed foundations, and pipe networks attached to structural elements. Understanding their vibratory behaviour is of interest for safer and more efficient design, in particular for better control of their response to moving loads or to stresses induced by external actions such as ambient vibrations, earthquakes, etc. Another field of application is the monitoring and control of the in-service condition of these structures using structural health monitoring (SHM) vibration methods.

Due to their wide range of applications, many studies have focused on dynamics of beams that rest on continuous foundations or on discrete flexible supports. For good reviews in this field one can refer to [1,2].

As for continuous foundations, the precursor works of [3–5] deal with Euler or Timoshenko beams, respectively, lying on an elastic or visco-elastic foundation of a Winkler compressional type. Let us mention the instability analysis developed by [6] and the work of [7,8] dealing with non-linear foundations and stratified foundations, respectively. The continuous approach has also been applied to discrete foundations by redistributing discrete forces in continuous form, e.g., [9].

Studies that explicitly take into account the discrete nature of supports were initiated by Mead for periodic [10] or non-periodic [11] distributions with the aim of determining their vibro-acoustic behaviour. More recently, ref. [2] investigates the steady-state dynamic
response under the influence of an axial force. It should be noted that in these studies the supports are often restricted to compression springs, and their rotational stiffness is much less frequently considered [10]. However, in a bending beam system, this type of spring necessarily plays an essential role that cannot be neglected in practice.

Several methods have been used to address the dynamic behaviour of multi-supported beams. Among them, the most usual are the transfer matrix method [9,12,13], the use of Green functions e.g., [2], and the numerical approach by the finite element method e.g., [7]. Asymptotic methods have been used less frequently. Particular mention should be made of the work of [14], who is interested in the dynamics of stretched beams with concentrated masses and discrete compressional elastic supports.

The objective of this paper is twofold: first, to develop a unified model that encompasses the phenomena of rotation/compression/bending within a synthetic formulation, thereby translating the overall behaviour of multi-supported structures; and second, to formulate the relevant boundary conditions in order to perform the modal analysis.

Despite recent advances in the field, to our knowledge, no unified model is available in the literature. Given the large number of possible situations, the study focuses on an Euler beam periodically supported by co-localised compression and rotation elastic springs, as illustrated in Figure 1. The stiffnesses of these springs, in comparison to that of the beam, can vary over a wide range. From the academic point of view, this example is both simple enough to allow analytical developments and rich enough to illustrate how to overcome some theoretical issues. Nevertheless, the framework of the method would allow Timoshenko beams or viscoelastic supports to be considered.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Above, irreducible period $\Omega$. Bottom, beam periodically supported by the spring system $(C, \mathcal{K})$.</figcaption></figure>

To establish this model, we will use an asymptotic method (homogenization of discrete periodic media; HPDM) whose advantages are as follows: (i) the determination, within a rigorous framework, of explicit macroscopic models that describe the long-wavelength dynamics on the basis of period parameters and their contrasts; (ii) the identification of local and global efforts; (iii) the delimitation of the domains of validity of each obtained model; (iv) the significant savings in computational resources (compared with a full calculation) thanks to the model reduction thus achieved; and (v) the physical understanding of the various possible situations by retaining only the dominant effects on a macroscopic scale.

The modal analysis in the long-wave range requires, in addition to the identification of the equivalent macroscopic model, the determination of the appropriate boundary conditions. Let us briefly comment on these two aspects.

The homogenization method is based on the assumption of a strong separation of scale between the macroscopic phenomenon and the microscopic structure. It has been first developed for continuous media [15–17], then extended to discrete systems, such as reticulated structures or continuous networks by [14,18–20]. The homogenization of periodic discrete media (HPDM) has in particular been used to establish generalised beam models [21–23]. Furthermore, homogenization analysis has also shown its interest in exploring
response under the influence of an axial force. It should be noted that in these studies the supports are often restricted to compression springs, and their rotational stiffness is much less frequently considered [10]. However, in a bending beam system, this type of spring necessarily plays an essential role that cannot be neglected in practice.

Several methods have been used to address the dynamic behaviour of multi-supported beams. Among them, the most usual are the transfer matrix method [9,12,13], the use of Green functions e.g., [2], and the numerical approach by the finite element method e.g., [7]. Asymptotic methods have been used less frequently. Particular mention should be made of the work of [14], who is interested in the dynamics of stretched beams with concentrated masses and discrete compressional elastic supports.

The objective of this paper is twofold: first, to develop a unified model that encompasses the phenomena of rotation/compression/bending within a synthetic formulation, thereby translating the overall behaviour of multi-supported structures; and second, to formulate the relevant boundary conditions in order to perform the modal analysis.

Despite recent advances in the field, to our knowledge, no unified model is available in the literature. Given the large number of possible situations, the study focuses on an Euler beam periodically supported by co-localised compression and rotation elastic springs, as illustrated in Figure 1. The stiffnesses of these springs, in comparison to that of the beam, can vary over a wide range. From the academic point of view, this example is both simple enough to allow analytical developments and rich enough to illustrate how to overcome some theoretical issues. Nevertheless, the framework of the method would allow Timoshenko beams or viscoelastic supports to be considered.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Above, irreducible period $\Omega$. Bottom, beam periodically supported by the spring system $(C, \mathcal{K})$.</figcaption></figure>

To establish this model, we will use an asymptotic method (homogenization of discrete periodic media; HPDM) whose advantages are as follows: (i) the determination, within a rigorous framework, of explicit macroscopic models that describe the long-wavelength dynamics on the basis of period parameters and their contrasts; (ii) the identification of local and global efforts; (iii) the delimitation of the domains of validity of each obtained model; (iv) the significant savings in computational resources (compared with a full calculation) thanks to the model reduction thus achieved; and (v) the physical understanding of the various possible situations by retaining only the dominant effects on a macroscopic scale.

The modal analysis in the long-wave range requires, in addition to the identification of the equivalent macroscopic model, the determination of the appropriate boundary conditions. Let us briefly comment on these two aspects.

The homogenization method is based on the assumption of a strong separation of scale between the macroscopic phenomenon and the microscopic structure. It has been first developed for continuous media [15–17], then extended to discrete systems, such as reticulated structures or continuous networks by [14,18–20]. The homogenization of periodic discrete media (HPDM) has in particular been used to establish generalised beam models [21–23]. Furthermore, homogenization analysis has also shown its interest in exploring
effects linked to imperfect scale separation through higher-order terms [24]. In quasi-static regimes, this work has made it possible to establish a link between phenomenological theories [25–27] of higher-order or micromorphic media and homogenised higher-order descriptions e.g., [28–30]. In the field of elastodynamics, these higher-order terms make it possible to tackle long wavelength dispersion, whether in 3D composite or stratified media, cf. [31–33].

The homogenised models are to be applied within the domain and must be completed at the edges of the domain. Nevertheless, the boundary conditions expressed on the variables of the homogenised model only imperfectly reflect the actual boundary conditions. This difficulty has been identified since the work of [34]. A way to overcome this drawback is to introduce a boundary layer that connects the true boundary condition to the homogenised condition inside the domain. The work of [35–37] provides a mathematical framework for these concepts. Despite these advances, this question is still topical, as shown by the recent work [38] in static regime, [39] in dynamic regime, and [40] for defining boundary layers in higher-order continuous media. As far as discrete media are concerned, to our knowledge, no specific work on this subject has been developed.

The study presented in this article deals with the various aspects mentioned above. The article is structured as follows: The first part is devoted to the determination of macroscopic models. These are constructed at the leading order using the HPDM method for each order of magnitude of stiffness contrasts in rotation and compression. The role of the contrasts on the order of the differential operator governing the overall dynamic behaviour is clearly highlighted. A so-called generic model, which encompasses all the descriptions, is established. It is shown by comparison with a Floquet-Bloch analysis that its range of validity extends well beyond the theoretical hypothesis of very large scale separation. The second part focuses on the formulation of relevant boundary conditions. A mixed "discrete/continuous" approach is proposed, which gives the effective boundary conditions to be applied to the macroscopic variables in order to reproduce the real conditions. The Robin-type conditions thus established associated with the macroscopic models can be used directly to carry out the modal analysis of these systems. The results obtained are very close to those obtained from a complete finite element calculation.

## 2. Method

In order to determine the eigenmodes of the periodic structure from the macroscopic model, the general procedure proposed is as follows:

1. Reduce the continuous structure studied to an equivalent discrete model. This is done by condensing the forces at the period nodes.

2. Derive the exact finite difference equation for the equilibrium between adjacent nodes in the discrete system.

3. Define the small scale parameter $\epsilon$, which is used (i) to define the orders of magnitude of the spring stiffnesses as a function of the bending stiffness of the beams and (ii) to develop the kinematic variables asymptotically.

4. Solve the problems in increasing order of $\epsilon$ until you get a macroscopic model where the order of the differential operator allows all physical boundary conditions to be considered.

5. A boundary condition relating to a kinematic variable is applied directly to the corresponding macroscopic variable. A force boundary condition, on the other hand, is applied by re-expressing the local equilibrium at the node under consideration using the macroscopic kinematic variables.

Note that the Floquet-Bloch theory provides the dispersion relation for a periodic material at any wavelength. However, knowledge of these properties is not sufficient to carry out a modal analysis of a domain of finite size, which requires the boundary conditions to be taken into account.
effects linked to imperfect scale separation through higher-order terms [24]. In quasi-static regimes, this work has made it possible to establish a link between phenomenological theories [25–27] of higher-order or micromorphic media and homogenised higher-order descriptions e.g., [28–30]. In the field of elastodynamics, these higher-order terms make it possible to tackle long wavelength dispersion, whether in 3D composite or stratified media, cf. [31–33].

The homogenised models are to be applied within the domain and must be completed at the edges of the domain. Nevertheless, the boundary conditions expressed on the variables of the homogenised model only imperfectly reflect the actual boundary conditions. This difficulty has been identified since the work of [34]. A way to overcome this drawback is to introduce a boundary layer that connects the true boundary condition to the homogenised condition inside the domain. The work of [35–37] provides a mathematical framework for these concepts. Despite these advances, this question is still topical, as shown by the recent work [38] in static regime, [39] in dynamic regime, and [40] for defining boundary layers in higher-order continuous media. As far as discrete media are concerned, to our knowledge, no specific work on this subject has been developed.

The study presented in this article deals with the various aspects mentioned above. The article is structured as follows: The first part is devoted to the determination of macroscopic models. These are constructed at the leading order using the HPDM method for each order of magnitude of stiffness contrasts in rotation and compression. The role of the contrasts on the order of the differential operator governing the overall dynamic behaviour is clearly highlighted. A so-called generic model, which encompasses all the descriptions, is established. It is shown by comparison with a Floquet-Bloch analysis that its range of validity extends well beyond the theoretical hypothesis of very large scale separation. The second part focuses on the formulation of relevant boundary conditions. A mixed "discrete/continuous" approach is proposed, which gives the effective boundary conditions to be applied to the macroscopic variables in order to reproduce the real conditions. The Robin-type conditions thus established associated with the macroscopic models can be used directly to carry out the modal analysis of these systems. The results obtained are very close to those obtained from a complete finite element calculation.

## 2. Method

In order to determine the eigenmodes of the periodic structure from the macroscopic model, the general procedure proposed is as follows:

1. Reduce the continuous structure studied to an equivalent discrete model. This is done by condensing the forces at the period nodes.

2. Derive the exact finite difference equation for the equilibrium between adjacent nodes in the discrete system.

3. Define the small scale parameter $\epsilon$, which is used (i) to define the orders of magnitude of the spring stiffnesses as a function of the bending stiffness of the beams and (ii) to develop the kinematic variables asymptotically.

4. Solve the problems in increasing order of $\epsilon$ until you get a macroscopic model where the order of the differential operator allows all physical boundary conditions to be considered.

5. A boundary condition relating to a kinematic variable is applied directly to the corresponding macroscopic variable. A force boundary condition, on the other hand, is applied by re-expressing the local equilibrium at the node under consideration using the macroscopic kinematic variables.

Note that the Floquet-Bloch theory provides the dispersion relation for a periodic material at any wavelength. However, knowledge of these properties is not sufficient to carry out a modal analysis of a domain of finite size, which requires the boundary conditions to be taken into account.
## 2.1. Euler Beam Dynamics

Let us consider a Euler beam (width $a$, length $l$, depth $h$, inertia $I = a^3h/12$, cross-section $A = ah$, Young modulus $E$, volumic mass $\rho$) in a harmonic regime at angular frequency $\omega$ that denotes ($v(s)$, $\theta(s)$) the transversal displacement and rotation field, where $s$ is the abscissa along the beam (the time dependency $e^{i\omega t}$ is omitted) and linked by the relationship:

$$ \theta(s) = \frac{dv(s)}{ds} \qquad (1) $$

The shear force $T_{\text{beam}}(s)$ and the momentum $M_{\text{beam}}(s)$ follow the local balances:

$$ \frac{dT_{\text{beam}}(s)}{ds} = \rho A \omega^2 v(s) \qquad (2a) $$

$$ \frac{dM_{\text{beam}}(s)}{ds} + T_{\text{beam}}(s) = 0 \qquad (2b) $$

The bending behaviour law for small deformations reads:

$$ M_{\text{beam}}(s) = -EI \frac{d\theta(s)}{ds} \qquad (3) $$

Introducing (1)–(3) in the derivative of (2b) yields:

$$ \frac{d^4 v(s)}{ds^4} - \left(\frac{\gamma}{\ell}\right)^4 v(s) = 0 \quad \text{with} \quad \gamma = \sqrt[4]{\frac{\rho A \omega^2}{EI}} \quad \ell = \frac{2\pi\ell_b}{\lambda_b} \qquad (4) $$

where $\gamma$ and $\lambda_b$ are, respectively, the dimensionless bending wavenumber and the bending wavelength. Denoting by $m$ and $K_f$ the mass and bending stiffness of the beam of length $l$, respectively, we define the reference angular frequency $\omega_b$ as

$$ m = \rho A \ell; \quad K_f = \frac{12EI}{\ell^3}; \quad \omega_b = \sqrt{\frac{EI}{\rho A \ell^4}} = \sqrt{\frac{K_f}{12m}}; \quad \gamma^4 = \left(\frac{\omega}{\omega_b}\right)^2 \qquad (5) $$

Integrating (4) with the displacement/rotations conditions ($v^B$, $\theta^B$) and ($v^E$, $\theta^E$) at the extremities $B$ and $E$ of the beam (see Figure 2) provides:

$$ v(s) = C_1 \cos(\gamma \frac{s}{\ell}) + C_2 \sin(\gamma \frac{s}{\ell}) + C_3 \cosh(\gamma \frac{s}{\ell}) + C_4 \sinh(\gamma \frac{s}{\ell}) \qquad (6) $$

with

$$ \begin{cases} 2C_1 = \left(1 + \frac{6}{\gamma^2}f_3(\gamma)\right)v^B - \frac{6}{\gamma^2}f_6(\gamma)v^E + \frac{4\ell}{\gamma^2}f_2(\gamma)\theta^B + \frac{2\ell}{\gamma^2}f_5(\gamma)\theta^E \\ 2C_2 = -\frac{12}{\gamma^3}f_1(\gamma)v^B + \frac{12}{\gamma^3}f_4(\gamma)v^E + \frac{\ell}{\gamma}\left(1 - \frac{6}{\gamma^2}f_3(\gamma)\right)\theta^B - \frac{6\ell}{\gamma^3}f_6(\gamma)\theta^E \\ 2C_3 = \left(1 - \frac{6}{\gamma^2}f_3(\gamma)\right)v^B + \frac{6}{\gamma^2}f_6(\gamma)v^E - \frac{4\ell}{\gamma^2}f_2(\gamma)\theta^B - \frac{2\ell}{\gamma^2}f_5(\gamma)\theta^E \\ 2C_4 = \frac{12}{\gamma^3}f_1(\gamma)v^B - \frac{12}{\gamma^3}f_4(\gamma)v^E + \frac{\ell}{\gamma}\left(1 + \frac{6}{\gamma^2}f_3(\gamma)\right)\theta^B + \frac{6\ell}{\gamma^3}f_6(\gamma)\theta^E \end{cases} $$
## 2.1. Euler Beam Dynamics

Let us consider a Euler beam (width $a$, length $l$, depth $h$, inertia $I = a^3h/12$, cross-section $A = ah$, Young modulus $E$, volumic mass $\rho$) in a harmonic regime at angular frequency $\omega$ that denotes ($v(s)$, $\theta(s)$) the transversal displacement and rotation field, where $s$ is the abscissa along the beam (the time dependency $e^{i\omega t}$ is omitted) and linked by the relationship:

$$ \theta(s) = \frac{dv(s)}{ds} \qquad (1) $$

The shear force $T_{\text{beam}}(s)$ and the momentum $M_{\text{beam}}(s)$ follow the local balances:

$$ \frac{dT_{\text{beam}}(s)}{ds} = \rho A \omega^2 v(s) \qquad (2a) $$

$$ \frac{dM_{\text{beam}}(s)}{ds} + T_{\text{beam}}(s) = 0 \qquad (2b) $$

The bending behaviour law for small deformations reads:

$$ M_{\text{beam}}(s) = -EI \frac{d\theta(s)}{ds} \qquad (3) $$

Introducing (1)–(3) in the derivative of (2b) yields:

$$ \frac{d^4 v(s)}{ds^4} - \left(\frac{\gamma}{\ell}\right)^4 v(s) = 0 \quad \text{with} \quad \gamma = \sqrt[4]{\frac{\rho A \omega^2}{EI}} \ell = \frac{2\pi\ell}{\lambda_b} \qquad (4) $$

where $\gamma$ and $\lambda_b$ are, respectively, the dimensionless bending wavenumber and the bending wavelength. Denoting by $m$ and $K_f$ the mass and bending stiffness of the beam of length $l$, respectively, we define the reference angular frequency $\omega_b$ as

$$ m = \rho A l; \quad K_f = \frac{12EI}{l^3}; \quad \omega_b = \sqrt{\frac{EI}{\rho A l^4}} = \sqrt{\frac{K_f}{12m}}; \quad \gamma^4 = \left(\frac{\omega}{\omega_b}\right)^2 \qquad (5) $$

Integrating (4) with the displacement/rotations conditions ($v^B$, $\theta^B$) and ($v^E$, $\theta^E$) at the extremities $B$ and $E$ of the beam (see Figure 2) provides:

$$ v(s) = C_1 \cos(\gamma \frac{s}{l}) + C_2 \sin(\gamma \frac{s}{l}) + C_3 \cosh(\gamma \frac{s}{l}) + C_4 \sinh(\gamma \frac{s}{l}) \qquad (6) $$

with

$$ \begin{cases} 2C_1 = \left(1 + \frac{6}{\gamma^2}f_3(\gamma)\right)v^B - \frac{6}{\gamma^2}f_6(\gamma)v^E + \frac{4\ell}{\gamma^2}f_2(\gamma)\theta^B + \frac{2\ell}{\gamma^2}f_5(\gamma)\theta^E \\ 2C_2 = -\frac{12}{\gamma^3}f_1(\gamma)v^B + \frac{12}{\gamma^3}f_4(\gamma)v^E + \frac{\ell}{\gamma}\left(1 - \frac{6}{\gamma^2}f_3(\gamma)\right)\theta^B - \frac{6\ell}{\gamma^3}f_6(\gamma)\theta^E \\ 2C_3 = \left(1 - \frac{6}{\gamma^2}f_3(\gamma)\right)v^B + \frac{6}{\gamma^2}f_6(\gamma)v^E - \frac{4\ell}{\gamma^2}f_2(\gamma)\theta^B - \frac{2\ell}{\gamma^2}f_5(\gamma)\theta^E \\ 2C_4 = \frac{12}{\gamma^3}f_1(\gamma)v^B - \frac{12}{\gamma^3}f_4(\gamma)v^E + \frac{\ell}{\gamma}\left(1 + \frac{6}{\gamma^2}f_3(\gamma)\right)\theta^B + \frac{6\ell}{\gamma^3}f_6(\gamma)\theta^E \end{cases} $$
and:

$$
\begin{align}
f_1(\gamma) &= \frac{\sin(\gamma) \cosh(\gamma) + \cos(\gamma) \sinh(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^3}{12} & f_4(\gamma) &= \frac{\sinh(\gamma) + \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^3}{12} \nonumber \\
f_2(\gamma) &= \frac{\sin(\gamma) \cosh(\gamma) - \cos(\gamma) \sinh(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma}{4} & f_5(\gamma) &= \frac{\sinh(\gamma) - \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma}{2} \tag{7} \\
f_3(\gamma) &= \frac{\sinh(\gamma) \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^2}{6} & f_6(\gamma) &= \frac{\cosh(\gamma) - \cos(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^2}{6} \nonumber
\end{align}
$$

The functions the functions {$f_i(\gamma)$}, $i = 1, \dots, 6$ have the following properties:

$$
\lim_{\gamma \to 0} f_i(\gamma) = 1 \quad ; \quad \forall \gamma, \forall i \in [4; 6], \quad f_i(\gamma) \neq 0 \qquad (8)
$$

Introducing (6) into (3) provides the expressions of the forces $\underline{F} = \begin{pmatrix} T_{\text{beam}} \\ M_{\text{beam}}/\ell \end{pmatrix}$ exerted by the beam at its extremities ($\underline{F}^B$ at the beginning; $\underline{F}^E$ at the end) as a function of the displacement $\underline{U} = \begin{pmatrix} v \\ \ell\theta \end{pmatrix}$ at its extremities ($\underline{U}^B$ at the beginning; $\underline{U}^E$ at the end). They are written in matrix form as:

$$
\begin{align}
\underline{F}^B(\underline{U}^B, \underline{U}^E) &= K_f(\underline{H}^{BB} \cdot \underline{U}^B + \underline{H}^{BE} \cdot \underline{U}^E) \notag \\
\underline{F}^E(\underline{U}^B, \underline{U}^E) &= K_f(\underline{H}^{EB} \cdot \underline{U}^B + \underline{H}^{EE} \cdot \underline{U}^E) \tag{9}
\end{align}
$$

where the following elasto-inertial matrices have been introduced:

$$
\begin{align*}
\underline{H}^{BB} &= \begin{bmatrix} f_1(\gamma) & \frac{1}{2}f_3(\gamma) \\ \frac{1}{2}f_3(\gamma) & \frac{1}{3}f_2(\gamma) \end{bmatrix} &
\underline{H}^{EE} &= \begin{bmatrix} -f_1(\gamma) & \frac{1}{2}f_3(\gamma) \\ \frac{1}{2}f_3(\gamma) & -\frac{1}{3}f_2(\gamma) \end{bmatrix} \\
\underline{H}^{BE} &= \begin{bmatrix} -f_4(\gamma) & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & \frac{1}{6}f_5(\gamma) \end{bmatrix} &
\underline{H}^{EB} &= -T\underline{H}^{BE} = \begin{bmatrix} f_4(\gamma) & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & -\frac{1}{6}f_5(\gamma) \end{bmatrix}
\end{align*}
$$

<figure><img src="image_5.png" /><figcaption>Figure 2. Kinematic variables and forces at the nodes defined in the local frame of each beam element.</figcaption></figure>

2.2. Dynamics of Multi-Supported Beam—Finite Difference Formulation

The period $\Omega$ of the system is composed by a beam of size $\ell$ and a system of springs ($K, C$) of positive stiffnesses. At each node $x_n = n\ell$ of the system, the springs exert a restoring force $F_n^K = -Kv_n$ and a restoring moment $M_n^C = -C\theta_n$. Hence, we define the stiffness matrix $\underline{R}$ and the restoring forces $\underline{f}_n$ as:

$$
\underline{f}_n = -K_f \underline{R} \cdot \underline{U}_n \quad ; \quad K_f \underline{R} = \begin{bmatrix} \mathcal{K} & 0 \\ 0 & C/\ell^2 \end{bmatrix} = K_f \begin{bmatrix} \beta & 0 \\ 0 & \delta \end{bmatrix}
$$

Expressing, at each node $n$ the balance $T_n$ of the transverse force, and the balance $M_n/\ell$ of
the $\ell$-reduced momentum yields

$$
\left( \begin{array}{c} T_n \\ M_n / \ell \end{array} \right) = \underline{F}^E(\underline{U}_{n-1}, \underline{U}_n) - \underline{F}^B(\underline{U}_n, \underline{U}_{n+1}) + \underline{f}_n = \underline{0} \qquad (10)
$$
and:

$$
\begin{align}
f_1(\gamma) &= \frac{\sin(\gamma) \cosh(\gamma) + \cos(\gamma) \sinh(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^3}{12} & f_4(\gamma) &= \frac{\sinh(\gamma) + \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^3}{12} \nonumber \\
f_2(\gamma) &= \frac{\sin(\gamma) \cosh(\gamma) - \cos(\gamma) \sinh(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma}{4} & f_5(\gamma) &= \frac{\sinh(\gamma) - \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma}{2} \tag{7} \\
f_3(\gamma) &= \frac{\sinh(\gamma) \sin(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^2}{6} & f_6(\gamma) &= \frac{\cosh(\gamma) - \cos(\gamma)}{1 - \cos(\gamma) \cosh(\gamma)} \frac{\gamma^2}{6} \nonumber
\end{align}
$$

The functions the functions {$f_i(\gamma)$}, $i = 1, \dots, 6$ have the following properties:

$$
\lim_{\gamma \to 0} f_i(\gamma) = 1 \quad ; \quad \forall \gamma, \forall i \in [4; 6], \quad f_i(\gamma) \neq 0 \qquad (8)
$$

Introducing (6) into (3) provides the expressions of the forces $\underline{F} = \begin{pmatrix} T_{\text{beam}} \\ M_{\text{beam}}/\ell \end{pmatrix}$ exerted by the beam at its extremities ($\underline{F}^B$ at the beginning; $\underline{F}^E$ at the end) as a function of the displacement $\underline{U} = \begin{pmatrix} v \\ \ell\theta \end{pmatrix}$ at its extremities ($\underline{U}^B$ at the beginning; $\underline{U}^E$ at the end). They are written in matrix form as:

$$
\begin{align}
\underline{F}^B (\underline{U}^B, \underline{U}^E) &= K_f (\underline{H}^{BB} \cdot \underline{U}^B + \underline{H}^{BE} \cdot \underline{U}^E) \notag \\
\underline{F}^E (\underline{U}^B, \underline{U}^E) &= K_f (\underline{H}^{EB} \cdot \underline{U}^B + \underline{H}^{EE} \cdot \underline{U}^E) \tag{9}
\end{align}
$$

where the following elasto-inertial matrices have been introduced:

$$
\begin{align*}
\underline{H}^{BB} &= \begin{bmatrix} f_1(\gamma) & \frac{1}{2}f_3(\gamma) \\ \frac{1}{2}f_3(\gamma) & \frac{1}{3}f_2(\gamma) \end{bmatrix} &
\underline{H}^{EE} &= \begin{bmatrix} -f_1(\gamma) & \frac{1}{2}f_3(\gamma) \\ \frac{1}{2}f_3(\gamma) & -\frac{1}{3}f_2(\gamma) \end{bmatrix} \\
\underline{H}^{BE} &= \begin{bmatrix} -f_4(\gamma) & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & \frac{1}{6}f_5(\gamma) \end{bmatrix} &
\underline{H}^{EB} &= -T\underline{H}^{BE} = \begin{bmatrix} f_4(\gamma) & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & -\frac{1}{6}f_5(\gamma) \end{bmatrix}
\end{align*}
$$

<figure><img src="image_6.png" /><figcaption>Figure 2. Kinematic variables and forces at the nodes defined in the local frame of each beam element.</figcaption></figure>

2.2. Dynamics of Multi-Supported Beam—Finite Difference Formulation

The period $\Omega$ of the system is composed by a beam of size $\ell$ and a system of springs ($K, C$) of positive stiffnesses. At each node $x_n = n\ell$ of the system, the springs exert a restoring force $F_n^K = -Kv_n$ and a restoring moment $M_n^C = -C\theta_n$. Hence, we define the stiffness matrix $\underline{R}$ and the restoring forces $\underline{f}_n$ as:

$$
\underline{f}_n = -K_f \underline{R} \cdot \underline{U}_n \quad ; \quad K_f \underline{R} = \begin{bmatrix} \mathcal{K} & 0 \\ 0 & C/\ell^2 \end{bmatrix} = K_f \begin{bmatrix} \beta & 0 \\ 0 & \delta \end{bmatrix}
$$

Expressing, at each node $n$ the balance $T_n$ of the transverse force, and the balance $M_n/\ell$ of
the $\ell$-reduced momentum yields

$$
\begin{pmatrix} T_n \\ M_n / \ell \end{pmatrix} = \underline{F}^E (\underline{U}_{n-1}, \underline{U}_n) - \underline{F}^B (\underline{U}_n, \underline{U}_{n+1}) + \underline{f}_n = \underline{0} \quad (10)
$$
which can be rewritten:

$$ \underline{\underline{H}}^{EB}.\underline{U}_{n-1} + (\underline{\underline{H}}^{EE} - \underline{\underline{H}}^{BB} - \underline{\underline{R}}).\underline{U}_{n} + \underline{T} \underline{\underline{H}}^{EB}.\underline{U}_{n+1} = \underline{0} \quad (11) $$

It is convenient to introduce the diagonal matrices $\underline{\underline{F}}(\omega) = \underline{\underline{H}}^{BB} - \underline{\underline{H}}^{EE}$ and $\underline{\underline{D}}(\omega) = \frac{1}{2}(\underline{\underline{H}}^{EB} + \underline{T} \underline{\underline{H}}^{EB})$ then the antisymmetric matrix $\underline{\underline{C}}(\omega) = \frac{1}{2}(\underline{\underline{H}}^{EB} - \underline{T} \underline{\underline{H}}^{EB})$, i.e., explicitly:

$$ \underline{\underline{D}}(\omega) = \begin{bmatrix} f_4(\gamma) & 0 \\ 0 & -\frac{1}{6}f_5(\gamma) \end{bmatrix}; \underline{\underline{C}}(\omega) = \begin{bmatrix} 0 & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & 0 \end{bmatrix}; \underline{\underline{F}}(\omega) = \begin{bmatrix} f_1(\gamma) & 0 \\ 0 & \frac{1}{3}f_2(\gamma) \end{bmatrix} $$

to rearrange (10) in a finite difference set. The two balance equations are coupled through the antisymmetric matrix $\underline{\underline{C}}(\omega)$:

$$ \frac{1}{K_f} \begin{pmatrix} T_n \\ M_n / \ell \end{pmatrix} = \underline{\underline{D}}(\omega) \cdot (\underline{\underline{U}}_{n-1} - 2\underline{\underline{U}}_n + \underline{\underline{U}}_{n+1}) + \underline{\underline{C}}(\omega) \cdot (\underline{\underline{U}}_{n-1} - \underline{\underline{U}}_{n+1}) + (2\underline{\underline{D}}(\omega) - 2\underline{\underline{F}}(\omega) - \underline{\underline{R}})\underline{\underline{U}}_n = \underline{\underline{0}} \quad (12) $$

## 2.3. Asymptotic Method

The switch from the finite difference description to a continuous description is made using the homogenization method for discrete media. This requires a net scale separation between the characteristic vibration length $L$ and the period length $\ell$, i.e., $\ell = \varepsilon L$, with $\varepsilon \ll 1$. Consequently, each local beam undergoes a quasi-static regime with a small dimensionless bending number $\gamma \ll 1$.

**Notation:** In the following, comparisons of orders of magnitude are made. The following notations are used:

* If an adimensional variable (e.g., $\gamma$) is of order $\varepsilon^m$, we note that

$$ \gamma = \gamma^* \varepsilon^m = O(\varepsilon^m) \text{ with } \gamma^* = O(1) $$

* an expression $A = O(\varepsilon) \Leftrightarrow A(1 + O(\varepsilon^1)) = 0$ means that all terms smaller than $\varepsilon$ with respect to $A$ are neglected. With this notation, it is not necessary to specify the dimensions of the terms to be compared.

The steps of the method are as follows:

**Continuation of discrete variables:** introduction of a continuous displacement $\underline{\underline{U}}(x) = \begin{pmatrix} v(x) \\ \ell\theta(x) \end{pmatrix}$ as in the nodes of the system:

$$ \underline{\underline{U}}(x_n = n\ell) = \underline{\underline{U}}_n = \begin{pmatrix} v_n \\ \ell\theta_n \end{pmatrix} \quad (13) $$

**Expansion of kinematic variables:** Under the assumption of scale separation, the variation in displacements between two consecutive nodes is small. It is then possible to expand the displacements at the neighbouring nodes of the $n^{th}$ node into Taylor series, giving rise to macroscopic derivatives denoted by $\partial_x$:

$$ \underline{\underline{U}}_{n \pm 1} = \underline{\underline{U}}(x_n \pm \varepsilon L) = \underline{\underline{U}}(x_n) \pm \varepsilon L \partial_x \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^2}{2} \partial_x^2 \underline{\underline{U}}(x_n) \pm \frac{(\varepsilon L)^3}{6} \partial_x^3 \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^4}{24} \partial_x^4 \underline{\underline{U}}(x_n) + O(\varepsilon^5) $$

Thus, the two finite difference terms of (12) become in the continuous representation:

$$ \underline{\underline{U}}_{n-1} - 2\underline{\underline{U}}_n + \underline{\underline{U}}_{n+1} = (\varepsilon L)^2 \partial_x^2 \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^4}{12} \partial_x^4 \underline{\underline{U}}(x_n) + O(\varepsilon^6) $$

$$ \underline{\underline{U}}_{n-1} - \underline{\underline{U}}_{n+1} = -2\varepsilon L \partial_x \underline{\underline{U}}(x_n) - \frac{(\varepsilon L)^3}{3} \partial_x^3 \underline{\underline{U}}(x_n) + O(\varepsilon^5) $$
which can be rewritten:

$$ \underline{\underline{H}}^{EB} \underline{U}_{n-1} + (\underline{\underline{H}}^{EE} - \underline{\underline{H}}^{BB} - \underline{\underline{R}}) \underline{U}_n + \underline{T} \underline{\underline{H}}^{EB} \underline{U}_{n+1} = \underline{0} \quad (11) $$

It is convenient to introduce the diagonal matrices $\underline{\underline{F}}(\omega) = \underline{\underline{H}}^{BB} - \underline{\underline{H}}^{EE}$ and $\underline{\underline{D}}(\omega) = \frac{1}{2}(\underline{\underline{H}}^{EB} + \underline{T} \underline{\underline{H}}^{EB})$ then the antisymmetric matrix $\underline{\underline{C}}(\omega) = \frac{1}{2}(\underline{\underline{H}}^{EB} - \underline{T} \underline{\underline{H}}^{EB})$, i.e., explicitly:

$$ \underline{\underline{D}}(\omega) = \begin{bmatrix} f_4(\gamma) & 0 \\ 0 & -\frac{1}{6}f_5(\gamma) \end{bmatrix}; \underline{\underline{C}}(\omega) = \begin{bmatrix} 0 & \frac{1}{2}f_6(\gamma) \\ -\frac{1}{2}f_6(\gamma) & 0 \end{bmatrix}; \underline{\underline{F}}(\omega) = \begin{bmatrix} f_1(\gamma) & 0 \\ 0 & \frac{1}{3}f_2(\gamma) \end{bmatrix} $$

to rearrange (10) in a finite difference set. The two balance equations are coupled through the antisymmetric matrix $\underline{\underline{C}}(\omega)$:

$$ \frac{1}{K_f} \begin{pmatrix} T_n \\ M_n / \ell \end{pmatrix} = \underline{\underline{D}}(\omega) \cdot (\underline{\underline{U}}_{n-1} - 2\underline{\underline{U}}_n + \underline{\underline{U}}_{n+1}) + \underline{\underline{C}}(\omega) \cdot (\underline{\underline{U}}_{n-1} - \underline{\underline{U}}_{n+1}) + (2\underline{\underline{D}}(\omega) - 2\underline{\underline{F}}(\omega) - \underline{\underline{R}})\underline{\underline{U}}_n = \underline{\underline{0}} \quad (12) $$

## 2.3. Asymptotic Method

The switch from the finite difference description to a continuous description is made using the homogenization method for discrete media. This requires a net scale separation between the characteristic vibration length $L$ and the period length $\ell$, i.e., $\ell = \varepsilon L$, with $\varepsilon \ll 1$. Consequently, each local beam undergoes a quasi-static regime with a small dimensionless bending number $\gamma \ll 1$.

**Notation:** In the following, comparisons of orders of magnitude are made. The following notations are used:

* If an adimensional variable (e.g., $\gamma$) is of order $\varepsilon^m$, we note that

$$ \gamma = \gamma^* \varepsilon^m = O(\varepsilon^m) \text{ with } \gamma^* = O(1) $$

* an expression $A = O(\varepsilon) \Leftrightarrow A(1 + O(\varepsilon^1)) = 0$ means that all terms smaller than $\varepsilon$ with respect to $A$ are neglected. With this notation, it is not necessary to specify the dimensions of the terms to be compared.

The steps of the method are as follows:

**Continuation of discrete variables:** introduction of a continuous displacement $\underline{\underline{U}}(x) = \begin{pmatrix} v(x) \\ \ell\theta(x) \end{pmatrix}$ as in the nodes of the system:

$$ \underline{\underline{U}}(x_n = n\ell) = \underline{\underline{U}}_n = \begin{pmatrix} v_n \\ \ell\theta_n \end{pmatrix} \quad (13) $$

**Expansion of kinematic variables:** Under the assumption of scale separation, the variation in displacements between two consecutive nodes is small. It is then possible to expand the displacements at the neighbouring nodes of the $n^{th}$ node into Taylor series, giving rise to macroscopic derivatives denoted by $\partial_x$:

$$ \underline{\underline{U}}_{n \pm 1} = \underline{\underline{U}}(x_n \pm \varepsilon L) = \underline{\underline{U}}(x_n) \pm \varepsilon L \partial_x \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^2}{2} \partial_x^2 \underline{\underline{U}}(x_n) \pm \frac{(\varepsilon L)^3}{6} \partial_x^3 \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^4}{24} \partial_x^4 \underline{\underline{U}}(x_n) + O(\varepsilon^5) $$

Thus, the two finite difference terms of (12) become in the continuous representation:

$$ \underline{\underline{U}}_{n-1} - 2\underline{\underline{U}}_n + \underline{\underline{U}}_{n+1} = (\varepsilon L)^2 \partial_x^2 \underline{\underline{U}}(x_n) + \frac{(\varepsilon L)^4}{12} \partial_x^4 \underline{\underline{U}}(x_n) + O(\varepsilon^6) $$

$$ \underline{\underline{U}}_{n-1} - \underline{\underline{U}}_{n+1} = -2\varepsilon L \partial_x \underline{\underline{U}}(x_n) - \frac{(\varepsilon L)^3}{3} \partial_x^3 \underline{\underline{U}}(x_n) + O(\varepsilon^5) $$
Moreover, the continuous kinematic variables are also asymptotically expanded into powers of $\epsilon$:

$$ \underline{U}(x) = \sum_{j \ge 0} \epsilon^j \underline{U}^{(j)}(x) \quad \text{with} \quad \underline{U}^{(j)}(x) = O(\epsilon^0 \underline{U}^{(0)}(x)) $$

**Expansion of frequencies:** as $\gamma \ll 1$, functions $\{f_i\}_{1 \le i \le 6}$ (see (7)) are expanded in a Taylor series:

$$ f_i(\gamma) = 1 + \sum_{j \ge 1} \gamma^j \partial_{\gamma}^j f_i(0) \quad (14) $$

Finally, the angular frequency is small compared to the reference frequency $\omega_b$ (see (5)), thus is expanded asymptotically around the zero pulsation:

$$ \omega = \gamma^2 \omega_b = 0 + \sum_{j \ge 1} \epsilon^j \omega^{(j)} \quad \text{with} \quad \omega^{(j)} = O(\epsilon^0 \omega_b) \quad (15) $$

By introducing (15) into (14), the matrices $\underline{D}, \underline{C}$ and $\underline{F}$ can be asymptotically expanded into powers of $\epsilon$.

**Normalization:** In this step, the geometric and mechanical contrasts are defined in order of magnitude of $\epsilon$. In accordance with the Euler-Bernoulli beam hypothesis, each beam element must be sufficiently slender (width $a \ll \ell$). This geometric contrast is expressed by considering (without loss of generality) that:

$$ a/\ell = e_0\epsilon = O(\epsilon) $$

In addition, the order of magnitude of the spring stiffnesses in translation and rotation is weighed in power of $\epsilon$ compared to the bending stiffness of a local beam:

$$ \beta = \frac{K}{K_f} = O(\epsilon^p) \quad ; \quad \delta = \frac{C}{K_f \ell^2} = O(\epsilon^q) $$

Playing with these two orders of magnitude can generate different behaviours in the medium on a macroscopic scale. This will be the subject of Section 3.

**Resolution:** All the previous expansions are introduced in (12), so that:

$$ \frac{1}{K_f} \left( \begin{array}{c} \mathbf{T}_n \\ \mathbf{M}_n / \ell \end{array} \right) = - \sum_{j \ge 0} \epsilon^j \left( \begin{array}{c} \mathbf{T}^{(j)} \\ \mathbf{M}^{(j)} / \ell \end{array} \right) = 0 \quad ; \quad \forall \epsilon \to 0 \quad (16) $$

Considering the equilibria for each order separately, the following successive equations arise for the first four orders:

$$ \begin{aligned} \left( \begin{array}{c} \mathbf{T}^{(-1)} \\ \mathbf{M}^{(-1)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} \mathbf{T}^{(0)} \\ \mathbf{M}^{(0)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(1)}(x) + \left( 2(\underline{D}^{(0)} - \underline{F}^{(0)}) - \underline{R}^{(0)} \right) \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} \mathbf{T}^{(1)} \\ \mathbf{M}^{(1)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(2)}(x) + \sum_{i=0}^{1} \left( 2(\underline{D}^{(i)} - \underline{F}^{(i)}) - \underline{R}^{(i)} \right) \underline{U}^{(1-i)}(x) - 2L\underline{C}^{(0)} \partial_x \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} \mathbf{T}^{(2)} \\ \mathbf{M}^{(2)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(3)}(x) + \sum_{i=0}^{2} \left( 2(\underline{D}^{(i)} - \underline{F}^{(i)}) - \underline{R}^{(i)} \right) \underline{U}^{(2-i)}(x) - 2L \sum_{i=0}^{1} \underline{C}^{(i)} \partial_x \underline{U}^{(1-i)}(x) + L^2 \underline{D}^{(0)} \partial_x^2 \underline{U}^{(0)}(x) = \underline{0} \\ \dots & \end{aligned} $$

The explicit formulation of these equations can be found in Appendix A. In order to encompass in a single framework the different possible levels of stiffness contrast expressed
Moreover, the continuous kinematic variables are also asymptotically expanded into powers of $\epsilon$:

$$ \underline{U}(x) = \sum_{j \ge 0} \epsilon^j \underline{U}^{(j)}(x) \quad \text{with} \quad \underline{U}^{(j)}(x) = O(\epsilon^0 \underline{U}^{(0)}(x)) $$

**Expansion of frequencies:** as $\gamma \ll 1$, functions $\{f_i\}_{1 \le i \le 6}$ (see (7)) are expanded in a Taylor series:

$$ f_i(\gamma) = 1 + \sum_{j \ge 1} \gamma^j \partial_{\gamma}^j f_i(0) \quad (14) $$

Finally, the angular frequency is small compared to the reference frequency $\omega_b$ (see (5)), thus is expanded asymptotically around the zero pulsation:

$$ \omega = \gamma^2 \omega_b = 0 + \sum_{j \ge 1} \epsilon^j \omega^{(j)} \quad \text{with} \quad \omega^{(j)} = O(\epsilon^0 \omega_b) \quad (15) $$

By introducing (15) into (14), the matrices $\underline{D}$, $\underline{C}$ and $\underline{F}$ can be asymptotically expanded into powers of $\epsilon$.

**Normalization:** In this step, the geometric and mechanical contrasts are defined in order of magnitude of $\epsilon$. In accordance with the Euler-Bernoulli beam hypothesis, each beam element must be sufficiently slender (width $a \ll \ell$). This geometric contrast is expressed by considering (without loss of generality) that:

$$ a/\ell = e_0\epsilon = O(\epsilon) $$

In addition, the order of magnitude of the spring stiffnesses in translation and rotation is weighed in power of $\epsilon$ compared to the bending stiffness of a local beam:

$$ \beta = \frac{K}{K_f} = O(\epsilon^p) \quad ; \quad \delta = \frac{C}{K_f \ell^2} = O(\epsilon^q) $$

Playing with these two orders of magnitude can generate different behaviours in the medium on a macroscopic scale. This will be the subject of Section 3.

**Resolution:** All the previous expansions are introduced in (12), so that:

$$ \frac{1}{K_f} \left( \frac{T_n}{M_n / \ell} \right) = - \sum_{j \ge 0} \epsilon^j \left( \frac{T^{(j)}}{M^{(j)} / \ell} \right) = 0 \quad ; \quad \forall \epsilon \to 0 \quad (16) $$

Considering the equilibria for each order separately, the following successive equations arise for the first four orders:

$$ \begin{aligned} \left( \begin{array}{c} T^{(-1)} \\ M^{(-1)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} T^{(0)} \\ M^{(0)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(1)}(x) + \left( 2(\underline{D}^{(0)} - \underline{F}^{(0)}) - \underline{R}^{(0)} \right) \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} T^{(1)} \\ M^{(1)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(2)}(x) + \sum_{i=0}^{1} \left( 2(\underline{D}^{(i)} - \underline{F}^{(i)}) - \underline{R}^{(i)} \right) \underline{U}^{(1-i)}(x) - 2L\underline{C}^{(0)} \partial_x \underline{U}^{(0)}(x) = \underline{0} \\ \left( \begin{array}{c} T^{(2)} \\ M^{(2)} \end{array} \right) &= -\underline{R}^{(-1)} \underline{U}^{(3)}(x) + \sum_{i=0}^{2} \left( 2(\underline{D}^{(i)} - \underline{F}^{(i)}) - \underline{R}^{(i)} \right) \underline{U}^{(2-i)}(x) - 2L \sum_{i=0}^{1} \underline{C}^{(i)} \partial_x \underline{U}^{(1-i)}(x) + L^2 \underline{D}^{(0)} \partial_x^2 \underline{U}^{(0)}(x) = \underline{0} \\ \dots & \end{aligned} $$

The explicit formulation of these equations can be found in Appendix A. In order to encompass in a single framework the different possible levels of stiffness contrast expressed
by $\beta$ and $\delta$, these parameters are also developed in the form of $\beta = \sum_{j=-1}^{4} \epsilon^j \beta^{(j)}$ and $\delta = \sum_{j=-1}^{4} \epsilon^j \delta^{(j)}$. Thus, a $\epsilon^p$ contrast for $\beta$ will result in $\beta^{(p)} \neq 0$ and $\beta^{(j \neq p)} = 0$ (and similarly for a $\epsilon^q$ contrast for $\delta$)

### 3. Multi-Supported Beam Models

The aim of this section is to establish the equivalent macroscopic models for multi-supported beams corresponding to different pairs of spring stiffnesses. First, two basic models are presented that show the crucial influence of stiffness contrasts. It is then shown that by combining the characteristics of these two models, it is possible to build a generic model that covers all possible contrast situations. Finally, the validity of these macroscopic descriptions is analysed by comparing the dispersion curves obtained for these models with those obtained using the classical Floquet-Bloch approach.

#### 3.1. On the Level of Beam/Support Contrasts

Obviously, the stiffness contrast between the beam and the supports plays a crucial role. Indeed, in the case of vanishing support stiffness, the multi-beam system tends to behave as a single long beam, while at very high support stiffness, the multi-beam system tends to behave like a chain of clamped (and therefore independent) beams. The situations of interest lie between these two extreme cases.

In order for the supports to play an effective role, it is necessary for the restoring forces to act as a source in the force and moment balance. With a view to an assessment, the localised restoring forces can be redistributed along the length $\ell$ of the beam. This gives the following relationships:

$$O(\frac{dT_{\text{beam}}}{ds}) \le O(\frac{K_v}{\ell}) \quad ; \quad O(\frac{dM_{\text{beam}}}{ds}) \le O(\frac{C\theta}{\ell})$$

In the scale separation assumption, assuming that both $v$ and $\theta$ vary with the macroscopic length, one derives from the beam constitutive laws that $O(\frac{dT_{\text{beam}}}{ds}) = O(EI\frac{v}{L^4}) = O(\frac{K_f\ell^3}{L^4}v)$ and $O(\frac{dM_{\text{beam}}}{ds}) = O(EI\frac{\theta}{L^2}) = O(\frac{K_f\ell^3}{L^2}\theta)$. These assessments reported in the previous relations provide the lower bound of the contrasts below that support the effect vanishes, i.e.,

$$O(K_f \frac{\ell^3}{L^4}) \le O(\frac{K}{\ell}) \implies O(\beta) \ge \epsilon^4 \quad ; \quad O(K_f \frac{\ell^3}{L^2}) \le O(\frac{C}{\ell}) \implies O(\delta) \ge \epsilon^2$$

On the other hand, the contrast $\beta$ is also limited by an upper bound. In fact, a strong support imposes that the rotation and the curvatures vary at the local scale so that $O(\theta) = O(v/\ell)$ and $O(\frac{d\theta}{ds}) = O(v/\ell^2)$, while, according to the scale separation assumption, the transverse force and moment vary at the large scale. This gives the following maximum estimate of $\frac{dT_{\text{beam}}}{ds}$ that the restoring force can balance

$$O(\frac{dT_{\text{beam}}}{ds}) = O(EI\frac{v}{\ell^2L^2}) = O(\frac{K_f\ell}{L^2}v)$$

Hence

$$O(K_f \frac{\ell}{L^2}) \ge O(\frac{K}{\ell}) \implies O(\beta) \le \epsilon^2$$

Thus, the contrasts of interest are in the ranges $\epsilon^4 \le O(\beta) \le \epsilon^2$ and $\epsilon^2 \le O(\delta)$.

#### 3.2. Basic Models

The two basic models are established by considering the cases where the compressional supports are as stiff as possible, i.e., $O(\beta) = \epsilon^2$ and the rotational supports are either stiff or soft.
by $\beta$ and $\delta$, these parameters are also developed in the form of $\beta = \sum_{j=-1}^{4} \epsilon^j \beta^{(j)}$ and $\delta = \sum_{j=-1}^{4} \epsilon^j \delta^{(j)}$. Thus, a $\epsilon^p$ contrast for $\beta$ will result in $\beta^{(p)} \neq 0$ and $\beta^{(j \neq p)} = 0$ (and similarly for a $\epsilon^q$ contrast for $\delta$)

### 3. Multi-Supported Beam Models

The aim of this section is to establish the equivalent macroscopic models for multi-supported beams corresponding to different pairs of spring stiffnesses. First, two basic models are presented that show the crucial influence of stiffness contrasts. It is then shown that by combining the characteristics of these two models, it is possible to build a generic model that covers all possible contrast situations. Finally, the validity of these macroscopic descriptions is analysed by comparing the dispersion curves obtained for these models with those obtained using the classical Floquet-Bloch approach.

#### 3.1. On the Level of Beam/Support Contrasts

Obviously, the stiffness contrast between the beam and the supports plays a crucial role. Indeed, in the case of vanishing support stiffness, the multi-beam system tends to behave as a single long beam, while at very high support stiffness, the multi-beam system tends to behave like a chain of clamped (and therefore independent) beams. The situations of interest lie between these two extreme cases.

In order for the supports to play an effective role, it is necessary for the restoring forces to act as a source in the force and moment balance. With a view to an assessment, the localised restoring forces can be redistributed along the length $\ell$ of the beam. This gives the following relationships:

$$O\left(\frac{dT_{\text{beam}}}{ds}\right) \le O\left(\frac{K_v}{\ell}\right) \quad ; \quad O\left(\frac{dM_{\text{beam}}}{ds}\right) \le O\left(\frac{C\theta}{\ell}\right)$$

In the scale separation assumption, assuming that both $v$ and $\theta$ vary with the macroscopic length, one derives from the beam constitutive laws that $O\left(\frac{dT_{\text{beam}}}{ds}\right) = O(EI\frac{v}{L^4}) = O\left(\frac{K_f\ell^3}{L^4}v\right)$ and $O\left(\frac{dM_{\text{beam}}}{ds}\right) = O(EI\frac{\theta}{L^2}) = O\left(\frac{K_f\ell^3}{L^2}\theta\right)$. These assessments reported in the previous relations provide the lower bound of the contrasts below that support the effect vanishes, i.e.,

$$O\left(K_f \frac{\ell^3}{L^4}\right) \le O\left(\frac{K}{\ell}\right) \implies O(\beta) \ge \epsilon^4 \quad ; \quad O\left(K_f \frac{\ell^3}{L^2}\right) \le O\left(\frac{C}{\ell}\right) \implies O(\delta) \ge \epsilon^2$$

On the other hand, the contrast $\beta$ is also limited by an upper bound. In fact, a strong support imposes that the rotation and the curvatures vary at the local scale so that $O(\theta) = O(v/\ell)$ and $O(\frac{d\theta}{ds}) = O(v/\ell^2)$, while, according to the scale separation assumption, the transverse force and moment vary at the large scale. This gives the following maximum estimate of $\frac{dT_{\text{beam}}}{ds}$ that the restoring force can balance

$$O\left(\frac{dT_{\text{beam}}}{ds}\right) = O\left(EI\frac{v}{\ell^2 L^2}\right) = O\left(\frac{K_f \ell}{L^2} v\right)$$

Hence

$$O\left(K_f \frac{\ell}{L^2}\right) \ge O\left(\frac{K}{\ell}\right) \implies O(\beta) \le \epsilon^2$$

Thus, the contrasts of interest are in the ranges $\epsilon^4 \le O(\beta) \le \epsilon^2$ and $\epsilon^2 \le O(\delta)$.

#### 3.2. Basic Models

The two basic models are established by considering the cases where the compressional supports are as stiff as possible, i.e., $O(\beta) = \epsilon^2$ and the rotational supports are either stiff or soft.
### 3.2.1. Multi-Supported Beam with Stiff Rotational Support

Let us consider multi-beam supported by stiff supports such that

$$ \beta = \frac{\mathcal{K}}{K_f} = \beta^{(2)} \varepsilon^2 \quad ; \quad \delta = \frac{\mathcal{C}}{K_f \ell^2} = \delta^{(0)} \varepsilon^0 $$

The macroscopic description at leading order is obtained by solving problems (A2) to (A4) in the Appendix A) in the growing order of powers of $\varepsilon$. Problems (A1) and (A2) are trivially verified with $\theta^{(-1)} = 0$. At the order $\varepsilon^1$, (A3) gives:

$$ \left\{ \begin{array}{l} \mathbb{T}^{(1)} = 0 \\ \mathbb{M}^{(1)} = (\delta^{(0)} + 1) L \theta^{(0)} - L \partial_x v^{(0)} = 0 \end{array} \right. $$

from which we can deduce the kinematic relationship to the dominant order:

$$ \theta^{(0)} = \frac{1}{(\delta^{(0)} + 1)} \partial_x v^{(0)} \qquad (17) $$

The equilibrium equation for the shear force at order 2 (A4) gives:

$$ \mathbb{T}^{(2)} = L^2 \partial_x \theta^{(0)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} - L^2 \partial_x^2 v^{(0)} = 0 $$

By introducing (17) into $\mathbb{T}^{(2)}$, we obtain the macroscopic description at leading order, with $\omega = \varepsilon\omega^{(1)} + o(\varepsilon)$:

$$ (\mathcal{K} - m\omega^2)v(x) - \mathcal{C}_{eq}\partial_x^2 v(x) = \mathcal{O}(\varepsilon) \qquad (18a) $$

$$ \theta(x) = \frac{C_{eq}}{C} \partial_x v(x) + \mathcal{O}(\varepsilon) \qquad (18b) $$

where $C_{eq}$ represents the effective rotational stiffness, which combines in series the stiffness $C$ and stiffness $K_f \ell^2$ of the beam, i.e.,

$$ \frac{1}{C_{eq}} = \frac{1}{C} + \frac{1}{K_f \ell^2} \qquad (19) $$

The rotational spring $C$ plays an essential role in the macroscopic balance equation (18a), which differs fundamentally from that of a bending beam. $C$ also affects the kinematic relationship (18b) between the rotation $\theta$ and the displacement $v$. In addition, the compression spring imposes a cut-off pulsation $\omega_c$ defined by:

$$ \omega_c = \sqrt{\frac{\mathcal{K}}{m}} \qquad (20) $$

below which the differential operator (18a) driving $v$ is not propagative.

The macroscopic description (18a) can be interpreted in terms of the balance of a macroscopic transverse force $-\mathcal{T}(x)$ linked to the displacement gradient by an apparent modulus $\mu = \frac{C_{eq}}{\ell}$.

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{\mathcal{K}}{\ell} v(x) \qquad (21a) $$

$$ \mathcal{T}(x) = -\frac{C_{eq}}{\ell} \partial_x v(x) \qquad (21b) $$

Hence this model, denoted MSB-Stiff in the sequel behaves macroscopically as a kind of shear-beam for which (i) the effects of the local bending and of the rotational spring
### 3.2.1. Multi-Supported Beam with Stiff Rotational Support

Let us consider multi-beam supported by stiff supports such that

$$ \beta = \frac{\mathcal{K}}{K_f} = \beta^{(2)} \varepsilon^2 \quad ; \quad \delta = \frac{\mathcal{C}}{K_f \ell^2} = \delta^{(0)} \varepsilon^0 $$

The macroscopic description at leading order is obtained by solving problems (A2) to (A4) in the Appendix A) in the growing order of powers of $\varepsilon$. Problems (A1) and (A2) are trivially verified with $\theta^{(-1)} = 0$. At the order $\varepsilon^1$, (A3) gives:

$$ \left\{ \begin{array}{l} \mathbb{T}^{(1)} = 0 \\ \mathbb{M}^{(1)} = (\delta^{(0)} + 1) L \theta^{(0)} - L \partial_x v^{(0)} = 0 \end{array} \right. $$

from which we can deduce the kinematic relationship to the dominant order:

$$ \theta^{(0)} = \frac{1}{(\delta^{(0)} + 1)} \partial_x v^{(0)} \quad (17) $$

The equilibrium equation for the shear force at order 2 (A4) gives:

$$ \mathbb{T}^{(2)} = L^2 \partial_x \theta^{(0)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} - L^2 \partial_x^2 v^{(0)} = 0 $$

By introducing (17) into $\mathbb{T}^{(2)}$, we obtain the macroscopic description at leading order, with $\omega = \varepsilon\omega^{(1)} + o(\varepsilon)$:

$$ (\mathcal{K} - m\omega^2)v(x) - \mathcal{C}_{eq}\partial_x^2 v(x) = \mathcal{O}(\varepsilon) \quad (18a) $$

$$ \theta(x) = \frac{C_{eq}}{C} \partial_x v(x) + \mathcal{O}(\varepsilon) \quad (18b) $$

where $C_{eq}$ represents the effective rotational stiffness, which combines in series the stiffness $C$ and stiffness $K_f \ell^2$ of the beam, i.e.,

$$ \frac{1}{C_{eq}} = \frac{1}{C} + \frac{1}{K_f \ell^2} \quad (19) $$

The rotational spring $C$ plays an essential role in the macroscopic balance equation (18a), which differs fundamentally from that of a bending beam. $C$ also affects the kinematic relationship (18b) between the rotation $\theta$ and the displacement $v$. In addition, the compression spring imposes a cut-off pulsation $\omega_c$ defined by:

$$ \omega_c = \sqrt{\frac{\mathcal{K}}{m}} \quad (20) $$

below which the differential operator (18a) driving $v$ is not propagative.

The macroscopic description (18a) can be interpreted in terms of the balance of a macroscopic transverse force $-\mathcal{T}(x)$ linked to the displacement gradient by an apparent modulus $\mu = \frac{C_{eq}}{\ell}$.

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{\mathcal{K}}{\ell} v(x) \quad (21a) $$

$$ \mathcal{T}(x) = -\frac{C_{eq}}{\ell} \partial_x v(x) \quad (21b) $$

Hence this model, denoted MSB-Stiff in the sequel behaves macroscopically as a kind of shear-beam for which (i) the effects of the local bending and of the rotational spring
are encapsulated in the effective shear modulus $\mu$ and (ii) the effect of the compressional spring results in an effective lineic force $-\frac{K}{l}v(x)$.

### 3.2.2. Multi-Supported Beam with Soft Rotational Support

This model is obtained by reducing the stiffness of the rotational spring (compared to the previous case), so that:

$$ \beta = \frac{K}{K_f} = \beta^{(2)}\epsilon^2 \quad ; \quad \delta = \frac{C}{K_f l^2} = \delta^{(2)}\epsilon^2 $$

As before, the problems (A1) to (A2) are trivially verified with $\theta^{(-1)} = 0$. For the next order, (A3) gives:

$$ \begin{cases} \mathbb{T}^{(1)} = 0 \\ \mathbb{M}^{(1)} = L\theta^{(0)} - L\partial_x v^{(0)} = 0 \end{cases} $$

which leads to the classical kinematic relationship of a Euler beam:

$$ \theta^{(0)} = \partial_x v^{(0)} $$

At order $\epsilon^2$, the equilibrium equations (A4) are:

$$ \begin{cases} \mathbb{T}^{(2)} = \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} = 0 \\ \mathbb{M}^{(2)} = L\theta^{(1)} - L\partial_x v^{(1)} = 0 \end{cases} $$

$\mathbb{T}^{(2)}$ gives the expression for the first angular frequency term $\omega^{(1)}$, and $\mathbb{M}^{(2)}$ gives the same kinematic relationship as in the previous order. Thus:

$$ \omega^{(1)} = \sqrt{\beta^{(2)}} \omega_b > 0 \quad ; \quad \theta^{(1)} = \partial_x v^{(1)} $$

At the order $\epsilon^3$, (A5) gives:

$$ \begin{cases} \mathbb{T}^{(3)} = -\frac{\rho L^2}{E e_0^2} 2\omega^{(1)}\omega^{(2)}v^{(0)} = 0 \\ \mathbb{M}^{(3)} = L\theta^{(2)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210Ee_0^2}(\omega^{(1)})^2 \right] L\theta^{(0)} - L\partial_x v^{(2)} - \frac{13\rho L^2}{210Ee_0^2}(\omega^{(1)})^2 L\partial_x v^{(0)} = 0 \end{cases} $$

As $\omega^{(1)} \neq 0$, $\mathbb{T}^{(3)}$ imposes $\omega^{(2)} = 0$, and $\mathbb{M}^{(3)}$ provides a second-order correction to the kinematic relationship:

$$ \theta^{(2)} = \partial_x v^{(2)} + \left( \frac{\beta^{(2)}}{15} - \delta^{(2)} \right) \theta^{(0)} $$

Finally, the equilibrium (A6) in shear force at order $\epsilon^4$ is written:

$$ \mathbb{T}^{(4)} = -L^2 \delta^{(2)} \partial_x^2 v^{(0)} + \frac{L^4}{12} \partial_x^4 v^{(0)} - \left[ \frac{(\beta^{(2)})^2}{60} + \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(3)} \right] v^{(0)} = 0 $$

Since $\omega^{(2)} = 0$, the angular frequency expands to: $\omega = \epsilon\omega^{(1)} + \epsilon^3\omega^{(3)} + o(\epsilon^3)$ and consequently $\epsilon^4 2\omega^{(1)}\omega^{(3)} = \omega^2 - (\epsilon\omega^{(1)})^2 + o(\epsilon^4)$. By introducing this expression into $\mathbb{T}^{(4)}$, as well as all those resulting from the solutions to the previous problems, we obtain the following macroscopic description:

$$ (-m\omega^2 + K_{eq})v(x) - C\partial_x^2 v(x) + EI\ell\partial_x^4 v(x) = O(\epsilon) \quad (22a) $$
are encapsulated in the effective shear modulus $\mu$ and (ii) the effect of the compressional spring results in an effective lineic force $-\frac{K}{l}v(x)$.

### 3.2.2. Multi-Supported Beam with Soft Rotational Support

This model is obtained by reducing the stiffness of the rotational spring (compared to the previous case), so that:

$$ \beta = \frac{K}{K_f} = \beta^{(2)}\epsilon^2 \quad ; \quad \delta = \frac{C}{K_f l^2} = \delta^{(2)}\epsilon^2 $$

As before, the problems (A1) to (A2) are trivially verified with $\theta^{(-1)} = 0$. For the next order, (A3) gives:

$$ \begin{cases} \mathbb{T}^{(1)} = 0 \\ \mathbb{M}^{(1)} = L\theta^{(0)} - L\partial_x v^{(0)} = 0 \end{cases} $$

which leads to the classical kinematic relationship of a Euler beam:

$$ \theta^{(0)} = \partial_x v^{(0)} $$

At order $\epsilon^2$, the equilibrium equations (A4) are:

$$ \begin{cases} \mathbb{T}^{(2)} = \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} = 0 \\ \mathbb{M}^{(2)} = L\theta^{(1)} - L\partial_x v^{(1)} = 0 \end{cases} $$

$\mathbb{T}^{(2)}$ gives the expression for the first angular frequency term $\omega^{(1)}$, and $\mathbb{M}^{(2)}$ gives the same kinematic relationship as in the previous order. Thus:

$$ \omega^{(1)} = \sqrt{\beta^{(2)}} \omega_b > 0 \quad ; \quad \theta^{(1)} = \partial_x v^{(1)} $$

At the order $\epsilon^3$, (A5) gives:

$$ \begin{cases} \mathbb{T}^{(3)} = -\frac{\rho L^2}{E e_0^2} 2\omega^{(1)}\omega^{(2)}v^{(0)} = 0 \\ \mathbb{M}^{(3)} = L\theta^{(2)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210Ee_0^2}(\omega^{(1)})^2 \right] L\theta^{(0)} - L\partial_x v^{(2)} - \frac{13\rho L^2}{210Ee_0^2}(\omega^{(1)})^2 L\partial_x v^{(0)} = 0 \end{cases} $$

As $\omega^{(1)} \neq 0$, $\mathbb{T}^{(3)}$ imposes $\omega^{(2)} = 0$, and $\mathbb{M}^{(3)}$ provides a second-order correction to the kinematic relationship:

$$ \theta^{(2)} = \partial_x v^{(2)} + \left( \frac{\beta^{(2)}}{15} - \delta^{(2)} \right) \theta^{(0)} $$

Finally, the equilibrium (A6) in shear force at order $\epsilon^4$ is written:

$$ \mathbb{T}^{(4)} = -L^2 \delta^{(2)} \partial_x^2 v^{(0)} + \frac{L^4}{12} \partial_x^4 v^{(0)} - \left[ \frac{(\beta^{(2)})^2}{60} + \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(3)} \right] v^{(0)} = 0 $$

Since $\omega^{(2)} = 0$, the angular frequency expands to: $\omega = \epsilon\omega^{(1)} + \epsilon^3\omega^{(3)} + o(\epsilon^3)$ and consequently $\epsilon^4 2\omega^{(1)}\omega^{(3)} = \omega^2 - (\epsilon\omega^{(1)})^2 + o(\epsilon^4)$. By introducing this expression into $\mathbb{T}^{(4)}$, as well as all those resulting from the solutions to the previous problems, we obtain the following macroscopic description:

$$ \left(-m\omega^2 + K_{eq}\right)v(x) - C\partial_x^2 v(x) + EI\ell\partial_x^4 v(x) = O(\epsilon) \quad (22a) $$
$$ \theta(x) = \partial_x v(x) + O(\epsilon) \qquad (22b) $$

$K_{eq}$ is an effective stiffness that takes into account the stiffness $60K_f$ due to the beam in bending and $K$ of the support in compression. Note that as $K = \beta K_f = O(\epsilon^2 K_f)$, $K_{eq}$ corresponds at this order to a series association between these two elements.

$$ K_{eq} = K\left(1 - \frac{K}{60K_f}\right) = \frac{1}{\frac{1}{K} + \frac{1}{60K_f}} + o(\epsilon^2) \qquad (23) $$

As before, the presence of $K_{eq}$ introduces a cutoff angular frequency below which the differential operator (22a) is not propagative.

$$ \omega_{c,eq} = \sqrt{\frac{K_{eq}}{m}} < \omega_c \qquad (24) $$

We note that in the model described by (22a) and (22b), the kinematic relationship of the elementary beam is transferred to the macroscopic scale as well as the bending term, which is translated by the same fourth-order differential operator. This leads to the introduction of the macroscopic transverse force $\mathcal{T}(x)$ and moment $\mathcal{M}(x)$, which allow the model to be rewritten in terms of the transverse force and moment balances as shown below.

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{K_{eq}}{\ell} v(x) \qquad (25a) $$

$$ \partial_x \mathcal{M}(x) + \mathcal{T}(x) = -\frac{C}{\ell} \theta(x) \qquad (25b) $$

$$ \mathcal{M}(x) = -EI\partial_x^2 v(x) \qquad (25c) $$

$$ \theta(x) = \partial_x v(x) \qquad (25d) $$

This formulation emphasises that the presence of the springs ($K, C$) simply results in forces of volume $-K_{eq}v(x)$ and a moment of volume $-C\partial_x v(x)$, without modifying the inherent behaviour of the beams.

### 3.3. Multi-Supported Generic Beam Model

The basic models make it possible to build a single generic model that incorporates the characteristics of each. To do this, we keep the equilibrium Equations (25a) and (25b), and we replace the kinematic relationship between the macroscopic variables (25d) by (18b). This leads to the following set, where ($K_{eq}, C_{eq}$) are the equivalent stiffnesses (19) and (23):

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{K_{eq}}{\ell} v(x) $$

$$ \partial_x \mathcal{M}(x) + \mathcal{T}(x) = -\frac{C}{\ell} \theta(x) $$

$$ \mathcal{M}(x) = -EI\partial_x^2 v(x) $$

$$ \theta(x) = \frac{C_{eq}}{C} \partial_x v(x) $$

Combining these equations leads to the leading-order macroscopic behaviour equation, where the $\omega_{c,eq}$ cut-off angular frequency is defined as (24):

$$ m(\omega^2 - \omega_{c,eq}^2)v(x) + C_{eq}\partial_x^2 v(x) - EI\ell\partial_x^4 v(x) = 0 \qquad (26) $$
$$ \theta(x) = \partial_x v(x) + \mathcal{O}(\varepsilon) \qquad (22b) $$

$\mathcal{K}_{eq}$ is an effective stiffness that takes into account the stiffness $60K_f$ due to the beam in bending and $\mathcal{K}$ of the support in compression. Note that as $\mathcal{K} = \beta K_f = \mathcal{O}(\varepsilon^2 K_f)$, $\mathcal{K}_{eq}$ corresponds at this order to a series association between these two elements.

$$ \mathcal{K}_{eq} = \mathcal{K}\left(1 - \frac{\mathcal{K}}{60K_f}\right) = \frac{1}{\frac{1}{\mathcal{K}} + \frac{1}{60K_f}} + o(\varepsilon^2) \qquad (23) $$

As before, the presence of $\mathcal{K}_{eq}$ introduces a cutoff angular frequency below which the differential operator (22a) is not propagative.

$$ \omega_{c,eq} = \sqrt{\frac{\mathcal{K}_{eq}}{m}} < \omega_c \qquad (24) $$

We note that in the model described by (22a) and (22b), the kinematic relationship of the elementary beam is transferred to the macroscopic scale as well as the bending term, which is translated by the same fourth-order differential operator. This leads to the introduction of the macroscopic transverse force $\mathcal{T}(x)$ and moment $\mathcal{M}(x)$, which allow the model to be rewritten in terms of the transverse force and moment balances as shown below.

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{\mathcal{K}_{eq}}{\ell} v(x) \qquad (25a) $$

$$ \partial_x \mathcal{M}(x) + \mathcal{T}(x) = -\frac{C}{\ell} \theta(x) \qquad (25b) $$

$$ \mathcal{M}(x) = -EI\partial_x^2 v(x) \qquad (25c) $$

$$ \theta(x) = \partial_x v(x) \qquad (25d) $$

This formulation emphasises that the presence of the springs $(\mathcal{K}, \mathcal{C})$ simply results in forces of volume $-\mathcal{K}_{eq}v(x)$ and a moment of volume $-\mathcal{C}\partial_x v(x)$, without modifying the inherent behaviour of the beams.

### 3.3. Multi-Supported Generic Beam Model

The basic models make it possible to build a single generic model that incorporates the characteristics of each. To do this, we keep the equilibrium Equations (25a) and (25b), and we replace the kinematic relationship between the macroscopic variables (25d) by (18b). This leads to the following set, where $(\mathcal{K}_{eq}, \mathcal{C}_{eq})$ are the equivalent stiffnesses (19) and (23):

$$ \partial_x \mathcal{T}(x) = \rho A \omega^2 v(x) - \frac{\mathcal{K}_{eq}}{\ell} v(x) $$

$$ \partial_x \mathcal{M}(x) + \mathcal{T}(x) = -\frac{C}{\ell} \theta(x) $$

$$ \mathcal{M}(x) = -EI\partial_x^2 v(x) $$

$$ \theta(x) = \frac{C_{eq}}{C} \partial_x v(x) $$

Combining these equations leads to the leading-order macroscopic behaviour equation, where the $\omega_{c,eq}$ cut-off angular frequency is defined as (24):

$$ m(\omega^2 - \omega_{c,eq}^2)v(x) + \mathcal{C}_{eq}\partial_x^2 v(x) - EI\ell\partial_x^4 v(x) = 0 \qquad (26) $$
### 3.4. Dimensional Analysis

A dimensional analysis of (26) divided by $m\omega_b^2 = K_f/12$ highlights the dominant effects according to the order of magnitude of the stiffnesses ($K = \beta K_f$ and $C = \delta K_f \ell^2$). To do this, the displacement is rewritten in the form

$$v = v^r V^* \qquad (27)$$

where $v^r$ is the reference displacement and $V^*$ the adimensional displacement, $V^* = O(\varepsilon^0)$. In addition, we perform the change of variables $\hat{x} = x/L$ (where $L = \ell/\varepsilon \gg \ell$ is related to the wavelength $\lambda$ by the relation $L = \lambda/2\pi$). By defining the contrasts ($\beta = \beta^* \varepsilon^p, \delta = \delta^* \varepsilon^q$) (26) becomes:

$$\underbrace{\left(\frac{\omega}{\omega_b}\right)^2 V^*}_{\mathcal{T}_\omega} - 12 \underbrace{\frac{\beta^* \varepsilon^p}{1 + \frac{\beta^*}{60} \varepsilon^p} V^*}_{\mathcal{T}_\beta} + 12 \underbrace{\frac{\delta^* \varepsilon^{q+2}}{1 + \delta^* \varepsilon^q} \partial_x^2 V^*}_{\mathcal{T}_\delta} - \underbrace{\varepsilon^4 \partial_x^4 V^*}_{\mathcal{T}_b} = 0 \qquad (28)$$

The terms $\mathcal{T}_\beta, \mathcal{T}_\delta, \mathcal{T}_b$ represent respectively the effect of the elastic restoring of the compression springs, the rotation springs, and the bending of the beam elements. The term $\mathcal{T}_\omega$ represents the force of inertia of the system represented adimensionally by $(\omega/\omega_b)^2 = \gamma^4 \ll 1$. The dynamic balance (28) requires the $\mathcal{T}_\omega - \mathcal{T}_\beta$ term linked to $V^*$ to be balanced by the elastic terms $\mathcal{T}_\delta$ and/or $\mathcal{T}_b$. This makes it possible to determine, according to the order of magnitude ($K, C$), the associated model, which retains only the dominant terms of the same order and the order of magnitude of the angular frequency.

#### 3.4.1. Family of Models

The generic model (26) degenerates into less rich models depending on the order of magnitude of the different terms. Three main types of degenerate models (29)–(31) occur in the case of strong compression springs ($\beta = O(\varepsilon^2)$). These expressions are general and, depending on the contrast, the ($C_{eq}, K_{eq}$) parameters themselves may degenerate and the cut-off angular frequency $\omega_{c,eq}$ may disappear.

**Multi-supported beam with stiff rotational support (MSB-Stiff) ($q \le 1$):**

In that case $\mathcal{T}_\beta = O(\mathcal{T}_\delta) = O(\mathcal{T}_\omega) \gg \mathcal{T}_b$, thus the action of the rotational springs and internal bending, included in $C_{eq}$, balances the inertia and elastic restoring forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) + C_{eq}\partial_x^2 v(x) = O(\varepsilon) \qquad (29)$$

**Multi-supported beam with medium rotational support (MSB-Medium) ($q = 2$):**

In that case, all the terms are of the same order of magnitude. The action of the rotational springs and the bending of the beam elements are differentiated and balance the inertia and elastic restoring forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) + C\partial_x^2 v(x) - EI\ell\partial_x^4 v(x) = O(\varepsilon) \qquad (30)$$

**Multi-supported beam with soft rotational support (MSB-Soft) ($q \ge 3$):**

In that case, $\mathcal{T}_\omega - \mathcal{T}_\beta = O(\mathcal{T}_b) \gg O(\mathcal{T}_\delta)$ This means that the action of the rotational springs disappears at the dominant order, and only the global bending balances the inertia forces of the restoring elastic forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) - EI\ell\partial_x^4 v(x) = O(\varepsilon) \qquad (31)$$

Table 1 shows all macroscopic descriptions depending on stiffness contrasts.
### 3.4. Dimensional Analysis

A dimensional analysis of (26) divided by $m\omega_b^2 = K_f/12$ highlights the dominant effects according to the order of magnitude of the stiffnesses ($K = \beta K_f$ and $C = \delta K_f \ell^2$). To do this, the displacement is rewritten in the form

$$v = v^r V^* \qquad (27)$$

where $v^r$ is the reference displacement and $V^*$ the adimensional displacement, $V^* = O(\varepsilon^0)$. In addition, we perform the change of variables $\hat{x} = x/L$ (where $L = \ell/\varepsilon \gg \ell$ is related to the wavelength $\lambda$ by the relation $L = \lambda/2\pi$). By defining the contrasts ($\beta = \beta^*\varepsilon^p, \delta = \delta^*\varepsilon^q$) (26) becomes:

$$\underbrace{\left(\frac{\omega}{\omega_b}\right)^2 V^*}_{\mathcal{T}_\omega} - 12 \underbrace{\frac{\beta^* \varepsilon^p}{1 + \frac{\beta^*}{60} \varepsilon^p} V^*}_{\mathcal{T}_\beta} + 12 \underbrace{\frac{\delta^* \varepsilon^{q+2}}{1 + \delta^* \varepsilon^q} \partial_x^2 V^*}_{\mathcal{T}_\delta} - \underbrace{\varepsilon^4 \partial_x^4 V^*}_{\mathcal{T}_b} = 0 \qquad (28)$$

The terms $\mathcal{T}_\beta, \mathcal{T}_\delta, \mathcal{T}_b$ represent respectively the effect of the elastic restoring of the compression springs, the rotation springs, and the bending of the beam elements. The term $\mathcal{T}_\omega$ represents the force of inertia of the system represented adimensionally by $(\omega/\omega_b)^2 = \gamma^4 \ll 1$. The dynamic balance (28) requires the $\mathcal{T}_\omega - \mathcal{T}_\beta$ term linked to $V^*$ to be balanced by the elastic terms $\mathcal{T}_\delta$ and/or $\mathcal{T}_b$. This makes it possible to determine, according to the order of magnitude ($K, C$), the associated model, which retains only the dominant terms of the same order and the order of magnitude of the angular frequency.

#### 3.4.1. Family of Models

The generic model (26) degenerates into less rich models depending on the order of magnitude of the different terms. Three main types of degenerate models (29)–(31) occur in the case of strong compression springs ($\beta = O(\varepsilon^2)$). These expressions are general and, depending on the contrast, the ($C_{eq}, K_{eq}$) parameters themselves may degenerate and the cut-off angular frequency $\omega_{c,eq}$ may disappear.

**Multi-supported beam with stiff rotational support (MSB-Stiff) ($q \le 1$):**

In that case $\mathcal{T}_\beta = O(\mathcal{T}_\delta) = O(\mathcal{T}_\omega) \gg \mathcal{T}_b$, thus the action of the rotational springs and internal bending, included in $C_{eq}$, balances the inertia and elastic restoring forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) + C_{eq}\partial_x^2 v(x) = O(\varepsilon) \qquad (29)$$

**Multi-supported beam with medium rotational support (MSB-Medium) ($q = 2$):**

In that case, all the terms are of the same order of magnitude. The action of the rotational springs and the bending of the beam elements are differentiated and balance the inertia and elastic restoring forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) + C\partial_x^2 v(x) - EI\ell\partial_x^4 v(x) = O(\varepsilon) \qquad (30)$$

**Multi-supported beam with soft rotational support (MSB-Soft) ($q \ge 3$):**

In that case, $\mathcal{T}_\omega - \mathcal{T}_\beta = O(\mathcal{T}_b) \gg O(\mathcal{T}_\delta)$ This means that the action of the rotational springs disappears at the dominant order, and only the global bending balances the inertia forces of the restoring elastic forces of the compression springs.

$$m(\omega^2 - \omega_{c,eq}^2)v(x) - EI\ell\partial_x^4 v(x) = O(\varepsilon) \qquad (31)$$

Table 1 shows all macroscopic descriptions depending on stiffness contrasts.
<table><caption>**Table 1.** Set of macroscopic descriptions of multi-supported beams as a function of the order of magnitude of the stiffness constraints ($\beta, \delta$). The first column shows the three models (i) MSB-Stiff (in blue), (ii) MSB-Medium (in orange), and (iii) MSB-Soft (in green).</caption><thead><tr><th></th><th>β<sup>(p≤1)</sup></th><th>β<sup>(2)</sup></th><th>β<sup>(3)</sup></th><th>β<sup>(4)</sup></th><th>β<sup>(p≥5)</sup></th></tr></thead><tbody><tr><td>δ<sup>(q≤-1)</sup></td><td rowspan="4">No scale separation</td><td>m(ω<sup>2</sup> - ω<sub>c</sub><sup>2</sup>)v(x) + K<sub>f</sub>ℓ<sup>2</sup>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + K<sub>f</sub>ℓ<sup>2</sup>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(0)</sup></td><td>m(ω<sup>2</sup> - ω<sub>c</sub><sup>2</sup>)v(x) + C<sub>eq</sub>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + C<sub>eq</sub>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(1)</sup></td><td>m(ω<sup>2</sup> - ω<sub>c</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + C∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(2)</sup></td><td>m(ω<sup>2</sup> - ω<sub>c,eq</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>m(ω<sup>2</sup> - ω<sub>c</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>mω<sup>2</sup>v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(q≥3)</sup></td><td>m(ω<sup>2</sup> - ω<sub>c,eq</sub><sup>2</sup>)v(x) - EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>m(ω<sup>2</sup> - ω<sub>c</sub><sup>2</sup>)v(x) - EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>mω<sup>2</sup>v(x) - EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td></td></tr></tbody></table>

Some of the results presented in [14] are recovered, in particular the role of compression springs, which reveal a cut-off frequency. Our study emphasises the essential role of rotation springs, which was not considered in [14].

### 3.4.2. Order of Magnitude of Frequencies

Denoting $\max(\mathcal{T}_\delta, \mathcal{T}_b) = O(\varepsilon^m)$, (with $m \ge 2$ by construction) there are three distinct situations for the order of magnitude of the frequency:

1. If $\mathcal{T}_\beta \ll \max(\mathcal{T}_\delta, \mathcal{T}_b)$, then $\mathcal{T}_\omega = O(\max(\mathcal{T}_\delta, \mathcal{T}_b))$. In these cases, the macroscopic descriptions are always propagative, and the angular frequency expands as follows:

$$ \omega^2 = O(\omega_b^2 \varepsilon^m) \qquad (32) $$

2. If $\mathcal{T}_\beta = O(\max(\mathcal{T}_\delta, \mathcal{T}_b))$, then $\mathcal{T}_\omega = O(\varepsilon^m)$ and the macroscopic behaviour equations are propagative only for angular frequencies (32) greater than $\omega_c$.

3. If $\mathcal{T}_\beta \gg \max(\mathcal{T}_\delta, \mathcal{T}_b)$, then it is necessary that the difference is $\mathcal{T}_\omega - \mathcal{T}_\beta = O(\varepsilon^m)$. In other words, the expansion of the angular frequency is of the form:

$$ \omega^2 - \omega_c^2 = O(\varepsilon^m \omega_b^2) \quad \text{where} \quad \left(\frac{\omega_c}{\omega_b}\right)^2 = O(\varepsilon^p) \qquad (33) $$

Macroscopic descriptions are only propagative for angular frequencies greater than a cut-off frequency $\omega_{c,eq}$.

### 3.5. Dispersion Relationships at Long Wavelengths

Macroscopic models assume that $\varepsilon \to 0$. To determine the extent of the validity domain of these models for finite values of $\varepsilon$, the dispersion relations of the various degenerate homogenised models (cf. Section 3.5.1) and of the generic model are compared with those obtained by the exact Floquet-Bloch analysis (cf. Section 3.5.2), which is not restricted by the scale separation condition.

#### 3.5.1. Dispersion Relationships for the Different Models

By inserting $v$ of the form $v(x) = V^*e^{ik_\omega x}$ in each of the three models (29)–(31) and (26), the following dispersion relationships are obtained:

**MSB-Stiff ($q \le 1$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta}(k_\omega \ell)^2 \qquad (34) $$

**MSB-Medium ($q = 2$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + 12\delta(k_\omega \ell)^2 + (k_\omega \ell)^4 \qquad (35) $$
<table><caption>**Table 1.** Set of macroscopic descriptions of multi-supported beams as a function of the order of magnitude of the stiffness constraints ($\beta, \delta$). The first column shows the three models (i) MSB-Stiff (in blue), (ii) MSB-Medium (in orange), and (iii) MSB-Soft (in green).</caption><thead><tr><th></th><th>β<sup>(p≤1)</sup></th><th>β<sup>(2)</sup></th><th>β<sup>(3)</sup></th><th>β<sup>(4)</sup></th><th>β<sup>(p≥5)</sup></th></tr></thead><tbody><tr><td>δ<sup>(q≤−1)</sup></td><td rowspan="4">No scale separation</td><td>m(ω<sup>2</sup>−ω<sub>c</sub><sup>2</sup>)v(x) + K<sub>f</sub>ℓ<sup>2</sup>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + K<sub>f</sub>ℓ<sup>2</sup>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(0)</sup></td><td>m(ω<sup>2</sup>−ω<sub>c</sub><sup>2</sup>)v(x) + C<sub>eq</sub>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + C<sub>eq</sub>∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(1)</sup></td><td>m(ω<sup>2</sup>−ω<sub>c</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td><td>mω<sup>2</sup>v(x) + C∂<sub>x</sub><sup>2</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(2)</sup></td><td>m(ω<sup>2</sup>−ω<sub>c,eq</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>m(ω<sup>2</sup>−ω<sub>c</sub><sup>2</sup>)v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>mω<sup>2</sup>v(x) + C∂<sub>x</sub><sup>2</sup>v(x) + EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td></td></tr><tr><td>δ<sup>(q≥3)</sup></td><td>m(ω<sup>2</sup>−ω<sub>c,eq</sub><sup>2</sup>)v(x) − EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>m(ω<sup>2</sup>−ω<sub>c</sub><sup>2</sup>)v(x) − EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td>mω<sup>2</sup>v(x) − EIℓ∂<sub>x</sub><sup>4</sup>v(x) = 0</td><td></td></tr></tbody></table>

Some of the results presented in [14] are recovered, in particular the role of compression springs, which reveal a cut-off frequency. Our study emphasises the essential role of rotation springs, which was not considered in [14].

### 3.4.2. Order of Magnitude of Frequencies

Denoting $\max(\mathcal{T}_\delta, \mathcal{T}_b) = O(\varepsilon^m)$, (with $m \ge 2$ by construction) there are three distinct situations for the order of magnitude of the frequency:

1. If $\mathcal{T}_\beta \ll \max(\mathcal{T}_\delta, \mathcal{T}_b)$, then $\mathcal{T}_\omega = O(\max(\mathcal{T}_\delta, \mathcal{T}_b))$. In these cases, the macroscopic descriptions are always propagative, and the angular frequency expands as follows:

$$ \omega^2 = O(\omega_b^2 \varepsilon^m) \qquad (32) $$

2. If $\mathcal{T}_\beta = O(\max(\mathcal{T}_\delta, \mathcal{T}_b))$, then $\mathcal{T}_\omega = O(\varepsilon^m)$ and the macroscopic behaviour equations are propagative only for angular frequencies (32) greater than $\omega_c$.

3. If $\mathcal{T}_\beta \gg \max(\mathcal{T}_\delta, \mathcal{T}_b)$, then it is necessary that the difference is $\mathcal{T}_\omega - \mathcal{T}_\beta = O(\varepsilon^m)$. In other words, the expansion of the angular frequency is of the form:

$$ \omega^2 - \omega_c^2 = O(\varepsilon^m \omega_b^2) \quad \text{where} \quad \left(\frac{\omega_c}{\omega_b}\right)^2 = O(\varepsilon^p) \qquad (33) $$

Macroscopic descriptions are only propagative for angular frequencies greater than a cut-off frequency $\omega_{c,eq}$.

### 3.5. Dispersion Relationships at Long Wavelengths

Macroscopic models assume that $\varepsilon \to 0$. To determine the extent of the validity domain of these models for finite values of $\varepsilon$, the dispersion relations of the various degenerate homogenised models (cf. Section 3.5.1) and of the generic model are compared with those obtained by the exact Floquet-Bloch analysis (cf. Section 3.5.2), which is not restricted by the scale separation condition.

#### 3.5.1. Dispersion Relationships for the Different Models

By inserting $v$ of the form $v(x) = V^*e^{ik_\omega x}$ in each of the three models (29)–(31) and (26), the following dispersion relationships are obtained:

**MSB-Stiff ($q \le 1$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta}(k_\omega \ell)^2 \qquad (34) $$

**MSB-Medium ($q = 2$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + 12\delta(k_\omega \ell)^2 + (k_\omega \ell)^4 \qquad (35) $$
**MSB-Soft ($q \ge 3$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + (k_\omega \ell)^4 \quad (36) $$

**MSB-Generic ($\forall q$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta}(k_\omega \ell)^2 + (k_\omega \ell)^4 \quad (37) $$

The wavenumber $k_\omega$ is real only above the cut-off angular frequency $\omega_{c,eq}$; otherwise, it is purely imaginary. These relationships are presented in Figure 3 for three multi-supported systems with $\beta$ and $\delta$ stiffness contrasts corresponding to each of the models, which correspond to a typical value of $\epsilon = 0.1$:

**Case ($\beta = 10^{-2}, \delta = 1$) (in blue):** model MSB-Stiff (34);

**Case ($\beta = 10^{-2}, \delta = 10^{-2}$) (in orange):** model MSB-Medium (35);

**Case ($\beta = 10^{-5}, \delta = 10^{-3}$) (in green):** model MSB-Stiff (36) in which $\omega_{c,eq} = 0$.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Dispersion relationships $\omega/\omega_b$ versus the adimensional wavenumber $k_\omega\ell/\pi$, $k_\omega \in [0, \pi/\ell]$. For the three types of contrast, the dispersion relationships derived from the models are plotted as solid lines. The relationships calculated by Floquet-Bloch for the same contrasts are in dots.</figcaption></figure>

Figure 3 also shows the dispersion curves obtained by the Floquet-Bloch analysis described below for the same contrast values.

### 3.5.2. Dispersion Relations by Floquet-Bloch

The exact dispersion relation for the system described by (12) is obtained by searching for $\underline{U}_n$ in the form:

$$ \underline{U}_n = \underline{U}_0 e^{-ik_\omega n\ell} \qquad k_\omega \in [0, \frac{\pi}{\ell}] \quad (38) $$

By introducing (38) into (12), we get:

$$ \underline{D} \underline{U}_0 = [2 \cos(k\ell) \underline{D}(\omega) + 2I \sin(k\ell) \underline{C}(\omega) - 2\underline{F}(\omega) - \underline{R}] \underline{U}_0 = \underline{0} $$
**MSB-Soft ($q \ge 3$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + (k_\omega \ell)^4 \quad (36) $$

**MSB-Generic ($\forall q$):**

$$ \left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta}(k_\omega \ell)^2 + (k_\omega \ell)^4 \quad (37) $$

The wavenumber $k_\omega$ is real only above the cut-off angular frequency $\omega_{c,eq}$; otherwise, it is purely imaginary. These relationships are presented in Figure 3 for three multi-supported systems with $\beta$ and $\delta$ stiffness contrasts corresponding to each of the models, which correspond to a typical value of $\varepsilon = 0.1$:

**Case ($\beta = 10^{-2}, \delta = 1$) (in blue):** model MSB-Stiff (34);

**Case ($\beta = 10^{-2}, \delta = 10^{-2}$) (in orange):** model MSB-Medium (35);

**Case ($\beta = 10^{-5}, \delta = 10^{-3}$) (in green):** model MSB-Stiff (36) in which $\omega_{c,eq} = 0$.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Dispersion relationships $\omega/\omega_b$ versus the adimensional wavenumber $k_\omega\ell/\pi$, $k_\omega \in [0, \pi/\ell]$. For the three types of contrast, the dispersion relationships derived from the models are plotted as solid lines. The relationships calculated by Floquet-Bloch for the same contrasts are in dots.</figcaption></figure>

Figure 3 also shows the dispersion curves obtained by the Floquet-Bloch analysis described below for the same contrast values.

### 3.5.2. Dispersion Relations by Floquet-Bloch

The exact dispersion relation for the system described by (12) is obtained by searching for $\underline{U}_n$ in the form:

$$ \underline{U}_n = \underline{U}_0 e^{-ik_\omega n\ell} \qquad k_\omega \in [0, \frac{\pi}{\ell}] \quad (38) $$

By introducing (38) into (12), we get:

$$ \underline{\underline{D}}.\underline{\underline{U}}_0 = [2 \cos(k\ell)\underline{\underline{D}}(\omega) + 2I \sin(k\ell)\underline{\underline{C}}(\omega) - 2\underline{\underline{F}}(\omega) - \underline{\underline{R}}] \underline{\underline{U}}_0 = \underline{\underline{0}} $$
The determinant of $\mathbb{D}$ must vanish so that $U_0$ is not identically zero. The dispersion relation is then written as $\det(\mathbb{D}) = 0$, i.e.:

$$ \mathcal{E}(\omega, k_{\omega}) = 0 = \frac{1}{3}(\beta + 2f_1(\gamma) - 2f_4(\gamma)\cos(k_{\omega}\ell))(3\delta + 2f_2(\gamma) + f_5(\gamma)\cos(k_{\omega}\ell)) - f_6(\gamma)^2\sin^2(k_{\omega}\ell) \quad (39) $$

This non-linear relationship is solved numerically for the contrasting stiffnesses $(\beta, \delta)$ indicated in Section 3.5.1. The results reported in Figure 3 show that the models MSB-Soft and MSB-Medium coincide with the exact calculations for wavenumbers $k_{\omega} \in [0, \frac{\pi}{\ell}]$. The MSB-Stiff model is tangent to the exact solution and is a good approximation of the latter up to $k_{\omega}\ell = \varepsilon \approx 0.3\pi$, i.e., wavelengths greater than three periods. These homogenised models at the dominant order are therefore valid well beyond the scale separation condition $\varepsilon \ll 1$.

This conclusion is supported by the use of the dispersion relation of the generic model. This is shown in Figure 4, where the dispersion from the generic model and that from the exact calculation are compared for the three types of contrast. For $k_{\omega} \in [0, \frac{\pi}{\ell}]$ the dispersion of the generic model coincides with that given by the exact calculations for the models MSB-Soft and MSB-Medium. Dispersion relationships are therefore not affected by model enrichment. Applying the generic model to the MSB-Stiff case makes it possible to follow not only the tangency but also the curvature of the exact solution, which leads to a good approximation up to $\varepsilon \approx 0.5\pi$.

In conclusion, the macroscopic models make it possible not only to recover (well beyond the limit of the scale separation hypothesis) the dispersion curves derived from the Floquet-Bloch approach but also to access the macroscopic forces (in a generalized sense), which are not accessible by Floquet-Bloch.

$$ (\beta = 10^{-2}; \delta = 1) \quad (\beta = 10^{-2}; \delta = 10^{-2}) \quad (\beta = 10^{-5}; \delta = 10^{-3}) $$

<figure><img src="image_7.png" /><figcaption>Figure 4. Dispersion relationships for $k_{\omega} \in [0, \frac{\pi}{\ell}]$ from the generic model (dashed line) for the three types of contrast, compared with the dispersions from the Floquet-Bloch calculation (in dots).</figcaption></figure>

## 4. Boundary Conditions

The macroscopic dynamic models established in the previous section have shown their validity in infinite media within the limit of wavelengths greater than three periods.

The modal analysis of such finite-length systems requires the specification of the boundary conditions. This question involves two difficulties:
The determinant of $\mathbb{D}$ must vanish so that $U_0$ is not identically zero. The dispersion relation is then written as $\det(\mathbb{D}) = 0$, i.e.:

$$ \mathcal{E}(\omega, k_\omega) = 0 = \frac{1}{3}(\beta + 2f_1(\gamma) - 2f_4(\gamma)\cos(k_\omega\ell))(3\delta + 2f_2(\gamma) + f_5(\gamma)\cos(k_\omega\ell)) - f_6(\gamma)^2\sin^2(k_\omega\ell) \quad (39) $$

This non-linear relationship is solved numerically for the contrasting stiffnesses $(\beta, \delta)$ indicated in Section 3.5.1. The results reported in Figure 3 show that the models MSB-Soft and MSB-Medium coincide with the exact calculations for wavenumbers $k_\omega \in [0, \frac{\pi}{\ell}]$. The MSB-Stiff model is tangent to the exact solution and is a good approximation of the latter up to $k_\omega\ell = \varepsilon \approx 0.3\pi$, i.e., wavelengths greater than three periods. These homogenised models at the dominant order are therefore valid well beyond the scale separation condition $\varepsilon \ll 1$.

This conclusion is supported by the use of the dispersion relation of the generic model. This is shown in Figure 4, where the dispersion from the generic model and that from the exact calculation are compared for the three types of contrast. For $k_\omega \in [0, \frac{\pi}{\ell}]$ the dispersion of the generic model coincides with that given by the exact calculations for the models MSB-Soft and MSB-Medium. Dispersion relationships are therefore not affected by model enrichment. Applying the generic model to the MSB-Stiff case makes it possible to follow not only the tangency but also the curvature of the exact solution, which leads to a good approximation up to $\varepsilon \approx 0.5\pi$.

In conclusion, the macroscopic models make it possible not only to recover (well beyond the limit of the scale separation hypothesis) the dispersion curves derived from the Floquet-Bloch approach but also to access the macroscopic forces (in a generalized sense), which are not accessible by Floquet-Bloch.

<figure><img src="image_8.png" /><figcaption>Figure 4. Dispersion relationships for $k_\omega \in [0, \frac{\pi}{\ell}]$ from the generic model (dashed line) for the three types of contrast, compared with the dispersions from the Floquet-Bloch calculation (in dots).</figcaption></figure>

## 4. Boundary Conditions

The macroscopic dynamic models established in the previous section have shown their validity in infinite media within the limit of wavelengths greater than three periods.

The modal analysis of such finite-length systems requires the specification of the boundary conditions. This question involves two difficulties:
* Is it possible (and if so, how?) to use the macroscopic force and/or displacement variables to express the boundary conditions actually imposed on the real system that condition the microscopic force and/or displacement variables? This is the subject of Section 4.1;

* If the order of the differential operator of the macroscopic model is less than that of the 4th order operator of the bending beams (4), the number of macroscopic conditions is less than that required to express the real local conditions. A higher-order model is constructed in Section 4.4 in order to obtain a well-posed problem that is compatible with the real boundary conditions.

### 4.1. Energy Approach

Let us consider a multi-supported beam consisting of N periods $\Omega$ (see Figure 1) with a total length of $H = N\ell$. The energy balance of the generic beam is established by multiplying (26) by $\frac{1}{2}v(x)$ and integrating over $[0, H]$. After two integrations by parts, we obtain:

$$ \frac{1}{2} \int_{0}^{H} m \omega^2 v(x)^2 dx = \frac{1}{2} \int_{0}^{H} \left( \ell \frac{M(x)^2}{EI} + \frac{(K_{\text{eq}}v(x))^2}{K_{\text{eq}}} + \frac{(C\theta(x))^2}{C_{\text{eq}}} \right) dx + W_{\text{BC}} + \mathcal{O}(\varepsilon) \quad (40) $$
$$ \text{where } W_{\text{BC}} = \frac{\ell}{2} \left[ \mathcal{T}(x)v(x) + M(x)\theta(x) \frac{C}{C_{\text{eq}}} \right]_{0}^{H} \quad \text{and} \quad \theta(x) = \frac{C_{\text{eq}}}{C} \partial_{x}v(x) $$

Equation (40) gives, at the leading order, the kinetic energy, the deformation energy, and the work $W_{BC}$ of the boundary conditions expressed in macroscopic variables.

### 4.2. Effective Boundary Conditions

#### 4.2.1. Displacement/Rotation Boundary Conditions

By construction in Section 2.3, the kinematic variables at the nodes are the macroscopic kinematic variables, see (13). Consequently, the kinematic boundary conditions imposed locally on the real system are directly related to the macroscopic kinematic variables. If at $x = 0$ (node 1) or $x = H$ (node $N+1$), the imposed kinematic (denoted $\underline{U}^{\text{imp}}$) conditions give the following macroscopic boundary conditions:

$$ \underline{U}(0) = \underline{U}_{1}^{\text{imp}} = \begin{pmatrix} v^{\text{imp}} \\ \ell\theta^{\text{imp}} \end{pmatrix}_{|x=0} \quad ; \quad \underline{U}(H) = \underline{U}_{N+1}^{\text{imp}} = \begin{pmatrix} v^{\text{imp}} \\ \ell\theta^{\text{imp}} \end{pmatrix}_{|x=H} \quad (41) $$

#### 4.2.2. Force/Moment Boundary Conditions

These forces, which are correctly described within the structure, are disturbed at the ends due to the absence of the adjacent node. In fact, at the extreme nodes, the equilibrium is no longer governed by (11) but by an expression that can be derived by (i) removing the contributions of elements outside the domain and adding the imposed forces; (ii) replacing the discrete variables at the nodes by continuous variables at these points according to (13). We therefore write at the extremities $x = 0$ and $x = H$:

$$ -\underline{H}^{\text{BB}}.\underline{U}(0) + \underline{H}^{\text{EB}}.\underline{U}(\ell) = \underline{F}^{\text{B}}\left(\begin{pmatrix} v(0) \\ \ell\theta(0) \end{pmatrix}, \begin{pmatrix} v(\ell) \\ \ell\theta(\ell) \end{pmatrix}\right) = \begin{pmatrix} T^{\text{imp}} \\ M^{\text{imp}}/\ell \end{pmatrix}_{|x=0} \quad (42) $$
$$ \underline{H}^{\text{EB}}.\underline{U}(H-\ell) + \left(\underline{H}^{\text{EE}} - \underline{R}\right).\underline{U}(H) = \underline{F}^{\text{E}}\left(\begin{pmatrix} v(H-\ell) \\ \ell\theta(H-\ell) \end{pmatrix}, \begin{pmatrix} v(H) \\ \ell\theta(H) \end{pmatrix}\right) - \underline{R} \cdot \begin{pmatrix} v(H) \\ \ell\theta(H) \end{pmatrix} = \begin{pmatrix} T^{\text{imp}} \\ M^{\text{imp}}/\ell \end{pmatrix}_{|x=H} $$

whose explicit formulation in terms of macroscopic kinematic variables is as follows at $x = 0$:

$$ v(0)f_1(\gamma) - v(\ell)f_4(\gamma) + \frac{\ell}{2}[\theta(0)f_3(\gamma) + \theta(\ell)f_6(\gamma)] = (T^{\text{imp}})_1 / K_f \quad (43a) $$

$$ \frac{1}{2}[v(0)f_3(\gamma) - v(\ell)f_6(\gamma)] + \frac{\ell}{3}[f_2(\gamma)\theta(0) + \frac{1}{2}\theta(\ell)f_5(\gamma)] = (M^{\text{imp}})_1 / (\ell K_f) \quad (43b) $$
* Is it possible (and if so, how?) to use the macroscopic force and/or displacement variables to express the boundary conditions actually imposed on the real system that condition the microscopic force and/or displacement variables? This is the subject of Section 4.1;

* If the order of the differential operator of the macroscopic model is less than that of the 4th order operator of the bending beams (4), the number of macroscopic conditions is less than that required to express the real local conditions. A higher-order model is constructed in Section 4.4 in order to obtain a well-posed problem that is compatible with the real boundary conditions.

### 4.1. Energy Approach

Let us consider a multi-supported beam consisting of N periods $\Omega$ (see Figure 1) with a total length of $H = N\ell$. The energy balance of the generic beam is established by multiplying (26) by $\frac{1}{2}v(x)$ and integrating over $[0, H]$. After two integrations by parts, we obtain:

$$ \frac{1}{2} \int_{0}^{H} m \omega^2 v(x)^2 dx = \frac{1}{2} \int_{0}^{H} \left( \ell \frac{M(x)^2}{EI} + \frac{(K_{\text{eq}}v(x))^2}{K_{\text{eq}}} + \frac{(C\theta(x))^2}{C_{\text{eq}}} \right) dx + W_{\text{BC}} + \mathcal{O}(\varepsilon) \quad (40) $$
$$ \text{where } W_{\text{BC}} = \frac{\ell}{2} \left[ \mathcal{T}(x)v(x) + M(x)\theta(x) \frac{C}{C_{\text{eq}}} \right]_{0}^{H} \quad \text{and} \quad \theta(x) = \frac{C_{\text{eq}}}{C} \partial_{x} v(x) $$

Equation (40) gives, at the leading order, the kinetic energy, the deformation energy, and the work $W_{BC}$ of the boundary conditions expressed in macroscopic variables.

### 4.2. Effective Boundary Conditions

#### 4.2.1. Displacement/Rotation Boundary Conditions

By construction in Section 2.3, the kinematic variables at the nodes are the macroscopic kinematic variables, see (13). Consequently, the kinematic boundary conditions imposed locally on the real system are directly related to the macroscopic kinematic variables. If at $x = 0$ (node 1) or $x = H$ (node $N+1$), the imposed kinematic (denoted $\underline{U}^{\text{imp}}$) conditions give the following macroscopic boundary conditions:

$$ \underline{U}(0) = \underline{U}_{1}^{\text{imp}} = \begin{pmatrix} v^{\text{imp}} \\ \ell\theta^{\text{imp}} \end{pmatrix}_{|x=0} \quad ; \quad \underline{U}(H) = \underline{U}_{N+1}^{\text{imp}} = \begin{pmatrix} v^{\text{imp}} \\ \ell\theta^{\text{imp}} \end{pmatrix}_{|x=H} \quad (41) $$

#### 4.2.2. Force/Moment Boundary Conditions

These forces, which are correctly described within the structure, are disturbed at the ends due to the absence of the adjacent node. In fact, at the extreme nodes, the equilibrium is no longer governed by (11) but by an expression that can be derived by (i) removing the contributions of elements outside the domain and adding the imposed forces; (ii) replacing the discrete variables at the nodes by continuous variables at these points according to (13). We therefore write at the extremities $x = 0$ and $x = H$:

$$ \begin{aligned} & -\underline{H}^{\text{BB}} \underline{U}(0) + \underline{H}^{\text{EB}} \underline{U}(\ell) = \underline{F}^{\text{B}} \left( \begin{pmatrix} v(0) \\ \ell\theta(0) \end{pmatrix}, \begin{pmatrix} v(\ell) \\ \ell\theta(\ell) \end{pmatrix} \right) = \begin{pmatrix} T^{\text{imp}} \\ M^{\text{imp}}/\ell \end{pmatrix}_{|x=0} \\ & \underline{H}^{\text{EB}} \underline{U}(H-\ell) + \left( \underline{H}^{\text{EE}} - \underline{R} \right) \underline{U}(H) = \underline{F}^{\text{E}} \left( \begin{pmatrix} v(H-\ell) \\ \ell\theta(H-\ell) \end{pmatrix}, \begin{pmatrix} v(H) \\ \ell\theta(H) \end{pmatrix} \right) - \underline{R} \cdot \begin{pmatrix} v(H) \\ \ell\theta(H) \end{pmatrix} = \begin{pmatrix} T^{\text{imp}} \\ M^{\text{imp}}/\ell \end{pmatrix}_{|x=H} \end{aligned} \quad (42) $$

whose explicit formulation in terms of macroscopic kinematic variables is as follows at $x = 0$:

$$ v(0)f_1(\gamma) - v(\ell)f_4(\gamma) + \frac{\ell}{2}[\theta(0)f_3(\gamma) + \theta(\ell)f_6(\gamma)] = (T^{\text{imp}})_1 / K_f \quad (43a) $$

$$ \frac{1}{2}[v(0)f_3(\gamma) - v(\ell)f_6(\gamma)] + \frac{\ell}{3}[f_2(\gamma)\theta(0) + \frac{1}{2}\theta(\ell)f_5(\gamma)] = (M^{\text{imp}})_1 / (\ell K_f) \quad (43b) $$
and at $x = H$

$$ -v(H)(f_1(\gamma) + \beta) + v(H - \ell)f_4(\gamma) + \frac{\ell}{2}[\theta(H)f_3(\gamma) + \theta(H - \ell)f_6(\gamma)] = (T^{\text{imp}})_{N+1}/K_f \quad (44a) $$

$$ \frac{1}{2}[v(H)f_3(\gamma) - v(H - \ell)f_6(\gamma)] - \frac{\ell}{3}[\theta(H)(f_2(\gamma) + \delta) + \frac{1}{2}\theta(H - \ell)f_5(\gamma)] = (M^{\text{imp}})_{N+1}/(K_f \ell) \quad (44b) $$

Note that this type of condition is similar to the Robin conditions coupled in $v$ and $\theta$. In particular, note that the nullity of the macroscopic forces $T(H)$ and $M(x)$ or of $W_{BC}$ does not reflect the free end condition. To put it in a nutshell:

1. A real condition on a kinematic variable is directly transposed into the homogenised model for the corresponding macroscopic kinematic variable.

2. A real condition on force or moment does not apply to the corresponding macroscopic effort of the homogenized model. In this case, it is necessary to express the force/moment equilibrium at the node using the macroscopic kinematic variables of this node and the adjacent internal node.

In all cases, this approach gives four boundary conditions expressed in macroscopic variables, which, for models of differential order four, make it possible to respect the conditions actually applied to a multi-supported beam of finite length. However, for models of differential order two, these four conditions are overabundant. This question is addressed in the next section.

### 4.3. Modal Analysis of a Differential Model of the Second Order

In this section, we focus on structures with stiff rotational springs that behave at the leading order as the model MSB-Stiff described by (29) (see Section 3.2.1). This situation needs to address the discrepancy between the number of boundary conditions required (two conditions at each edge beam) and the order of the differential operator in the macroscopic description (one condition at each extremity). In this view, three strategies are studied:

* Use the macroscopic model (29) whose differential operator is of order two and where only two boundary conditions can be expressed. The driving variable is the displacement, while the rotation is merely a hidden variable derived from it. This forces us to consider the boundary condition in displacement $v^{\text{imp}}$ and forces $T^{\text{imp}}$ (and disregard rotations $\theta^{\text{imp}}$ and moments $M^{\text{imp}}$) in (41) and (42).

* Use the extension of the model by the generic model (26) with a differential operator of the fourth order. In this case, four conditions (41) and (42) can be effectively applied.

* Use the model (55) developed in Section 4.4.2, which is enriched at higher orders and results in a fourth-order differential operator that allows the application of four boundary conditions (41) and (42).

These options have been explored in the context of the modal analysis of a multi-supported beam of kind MSB-Stiff ($\beta = 10^{-2}; \delta = 1$) with fixed-fixed and fixed-free boundary conditions.

#### 4.3.1. "Reduced" Macroscopic Modal Analysis

As an example, let us consider a fixed-fixed MSB-Stiff beam consisting of $N = 10$ periods. The "reduced" macroscopic modal analysis is based on (i) the MSB-Stiff model and (ii) the two zero displacement conditions imposed at the edges (and therefore ignores the zero rotation conditions imposed). This situation leads directly to the macroscopic modal deformation of the $j$th mode:

$$ v(x) = A \sin((k_{\omega})_j x); \quad \theta(x) = \frac{C_{\text{eq}}}{C} \partial_x v(x) = A(k_{\omega})_j \cos((k_{\omega})_j x), \quad \text{with } (k_{\omega})_j = \frac{j\pi}{H}, \quad j \in \mathbb{N}^* \quad (45) $$

where the angular frequency $\omega_j$, which is associated with $(k_\omega)_j$ is defined by the dispersion relation (34). Figure 5 shows the modal deformations of the first three modes thus obtained. For comparison, the numerical results from a finite element calculation (taking into account
and at $x = H$

$$ -v(H)(f_1(\gamma) + \beta) + v(H - \ell)f_4(\gamma) + \frac{\ell}{2}[\theta(H)f_3(\gamma) + \theta(H - \ell)f_6(\gamma)] = (T^{\text{imp}})_{N+1}/K_f \quad (44a) $$

$$ \frac{1}{2}[v(H)f_3(\gamma) - v(H - \ell)f_6(\gamma)] - \frac{\ell}{3}\left[\theta(H)(f_2(\gamma) + \delta) + \frac{1}{2}\theta(H - \ell)f_5(\gamma)\right] = (M^{\text{imp}})_{N+1}/(K_f \ell) \quad (44b) $$

Note that this type of condition is similar to the Robin conditions coupled in $v$ and $\theta$. In particular, note that the nullity of the macroscopic forces $T(H)$ and $M(x)$ or of $W_{BC}$ does not reflect the free end condition. To put it in a nutshell:

1. A real condition on a kinematic variable is directly transposed into the homogenised model for the corresponding macroscopic kinematic variable.

2. A real condition on force or moment does not apply to the corresponding macroscopic effort of the homogenized model. In this case, it is necessary to express the force/moment equilibrium at the node using the macroscopic kinematic variables of this node and the adjacent internal node.

In all cases, this approach gives four boundary conditions expressed in macroscopic variables, which, for models of differential order four, make it possible to respect the conditions actually applied to a multi-supported beam of finite length. However, for models of differential order two, these four conditions are overabundant. This question is addressed in the next section.

### 4.3. Modal Analysis of a Differential Model of the Second Order

In this section, we focus on structures with stiff rotational springs that behave at the leading order as the model MSB-Stiff described by (29) (see Section 3.2.1). This situation needs to address the discrepancy between the number of boundary conditions required (two conditions at each edge beam) and the order of the differential operator in the macroscopic description (one condition at each extremity). In this view, three strategies are studied:

* Use the macroscopic model (29) whose differential operator is of order two and where only two boundary conditions can be expressed. The driving variable is the displacement, while the rotation is merely a hidden variable derived from it. This forces us to consider the boundary condition in displacement $v^{\text{imp}}$ and forces $T^{\text{imp}}$ (and disregard rotations $\theta^{\text{imp}}$ and moments $M^{\text{imp}}$) in (41) and (42).

* Use the extension of the model by the generic model (26) with a differential operator of the fourth order. In this case, four conditions (41) and (42) can be effectively applied.

* Use the model (55) developed in Section 4.4.2, which is enriched at higher orders and results in a fourth-order differential operator that allows the application of four boundary conditions (41) and (42).

These options have been explored in the context of the modal analysis of a multi-supported beam of kind MSB-Stiff ($\beta = 10^{-2}$; $\delta = 1$) with fixed-fixed and fixed-free boundary conditions.

#### 4.3.1. "Reduced" Macroscopic Modal Analysis

As an example, let us consider a fixed-fixed MSB-Stiff beam consisting of $N = 10$ periods. The "reduced" macroscopic modal analysis is based on (i) the MSB-Stiff model and (ii) the two zero displacement conditions imposed at the edges (and therefore ignores the zero rotation conditions imposed). This situation leads directly to the macroscopic modal deformation of the $j$th mode:

$$ v(x) = A \sin((k_{\omega})_j x); \quad \theta(x) = \frac{C_{\text{eq}}}{C} \partial_x v(x) = A(k_{\omega})_j \cos((k_{\omega})_j x), \quad \text{with } (k_{\omega})_j = \frac{j\pi}{H}, \quad j \in \mathbb{N}^* \quad (45) $$

where the angular frequency $\omega_j$, which is associated with $(k_\omega)_j$ is defined by the dispersion relation (34). Figure 5 shows the modal deformations of the first three modes thus obtained. For comparison, the numerical results from a finite element calculation (taking into account
the zero rotation conditions) of the real fixed-fixed structure are also shown (after mesh study, eight bending beam elements are considered per beam).

In the model MSB-Stiff, the stiffness of the rotation springs is of the same order of magnitude as the bending stiffness of each beam element. The modal deformation of the system can be decomposed into (i) the amplitude of the displacements at the nodes and (ii) the deformation of the beams between the nodes. It can be seen that in the case of MSB-Stiff beams, the displacements (and rotations) of the nodes are fixed by the total length of the system, regardless of the stiffness of the supports and beams. On the other hand, the rotational supports play a significant role in the inter-node deformation of the beams. This effect is clearly visible both for rotational deformations (angular points at the nodes of the system caused by the discontinuity of moments between two elements) and for displacement deformations (curvature of each beam element linked to the rotation differential between two nodes). The compression springs do not influence the modal deformation, but they do influence the modal frequency; see (34).

There is good general agreement between the macroscopic model and the numerical model. The omission of the zero rotation condition in the model MSB-Stiff is clearly visible near the edges, mainly in the first period on each side. It should be noted that the differences between the homogenised model and the calculation increase from the fundamental mode to the third mode. In the latter case, however, they remain relatively limited due to the small-scale separation.

**Remark 1.** The choice of zero rotation would have led to $\theta(x) = A \sin((k_\omega)_j x)$ and therefore a deformation $v(x) = B_1 \cos((k_\omega)_j x) + B_2$. This expression is not admissible because it does not respect the fixed-fixed conditions on the driving variable. This is confirmed by the modal deformations obtained from the numerical simulations (see Figure 5), which show that $v$ is indeed of the form (45).

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

**Figure 5.** First three modal deformations of a 10-period fixed-fixed structure ($\beta = 10^{-2}; \delta = 1$) of the MSB-Stiff model (in blue) compared with the finite element model (in black). From left to right: (a) modal deformation in displacement $v(x)$, (b) difference between homogenised model and finite element model $v_{Macro}(x) - v_{FE}$, (c) modal deformation in rotation $\theta(x)$, (d) difference between homogenised model, and finite element model $\theta_{Macro}(x) - \theta_{FE}(x)$. The dots represent the nodes between the beams where the springs are located.
the zero rotation conditions) of the real fixed-fixed structure are also shown (after mesh study, eight bending beam elements are considered per beam).

In the model MSB-Stiff, the stiffness of the rotation springs is of the same order of magnitude as the bending stiffness of each beam element. The modal deformation of the system can be decomposed into (i) the amplitude of the displacements at the nodes and (ii) the deformation of the beams between the nodes. It can be seen that in the case of MSB-Stiff beams, the displacements (and rotations) of the nodes are fixed by the total length of the system, regardless of the stiffness of the supports and beams. On the other hand, the rotational supports play a significant role in the inter-node deformation of the beams. This effect is clearly visible both for rotational deformations (angular points at the nodes of the system caused by the discontinuity of moments between two elements) and for displacement deformations (curvature of each beam element linked to the rotation differential between two nodes). The compression springs do not influence the modal deformation, but they do influence the modal frequency; see (34).

There is good general agreement between the macroscopic model and the numerical model. The omission of the zero rotation condition in the model MSB-Stiff is clearly visible near the edges, mainly in the first period on each side. It should be noted that the differences between the homogenised model and the calculation increase from the fundamental mode to the third mode. In the latter case, however, they remain relatively limited due to the small-scale separation.

**Remark 1.** The choice of zero rotation would have led to $\theta(x) = A \sin((k_\omega)_j x)$ and therefore a deformation $v(x) = B_1 \cos((k_\omega)_j x) + B_2$. This expression is not admissible because it does not respect the fixed-fixed conditions on the driving variable. This is confirmed by the modal deformations obtained from the numerical simulations (see Figure 5), which show that $v$ is indeed of the form (45).

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

**Figure 5.** First three modal deformations of a 10-period fixed-fixed structure ($\beta = 10^{-2}; \delta = 1$) of the MSB-Stiff model (in blue) compared with the finite element model (in black). From left to right: (a) modal deformation in displacement $v(x)$, (b) difference between homogenised model and finite element model $v_{Macro}(x) - v_{FE}$, (c) modal deformation in rotation $\theta(x)$, (d) difference between homogenised model, and finite element model $\theta_{Macro}(x) - \theta_{FE}(x)$. The dots represent the nodes between the beams where the springs are located.
### 4.3.2. "Full" Macroscopic Modal Analysis

The full macroscopic modal analysis was carried out on a structure consisting of $N = 30$ periods. The two models used are that of the generic beam and that of the beam MSB-Stiff enriched with the higher orders; see Section 4.4. In each case, the four boundary conditions are those indicated in Section 4.2.

Figure 6 shows the fundamental modes determined by these two models with boundary conditions either fixed-fixed (Figure 6-top) or fixed-free (Figure 6-bottom). For comparison, the results obtained using finite elements and those obtained using "reduced" modal analysis are also shown.

Compared to the 10-period fixed-fixed structure, we can see that tripling the number of periods, as assumed, results in better scale separation and a gain in accuracy of the reduced modal analysis approach. Moreover, the concentration of the error for rotations on the periods at the edges is even more obvious.

For both types of boundary conditions, it can be seen that "full" modal analyses performed with four boundary conditions significantly improve the results. Furthermore, it appears that the higher-order model leads to a better approximation than the generic model. This is explained by the fact that the additional term in the higher-order model is inherited from effects specific to the contrast under study. This improves the description as it provides a model of an accuracy up to $O(\varepsilon^3)$. In contrast, the generic model is built up in such a way that it gives an accuracy up to $O(\varepsilon)$ regardless of the contrasts. It is always valid at the dominant order in the sense that the terms naturally vanish when they are not relevant. Nevertheless, even if the differential operator of the generic model is of fourth order, one cannot expect that it will reach the same level of accuracy as the higher-order model.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** First modal deformation of periodic structures ($N = 30$ periods, ($\beta = 10^{-2}$; $\delta = 1$)) with fixed-fixed (top) and fixed-free (bottom) boundary conditions. The homogenised models of order $O(\varepsilon)$ are shown in blue (MSB-Stiff) and green (MSB-Generic); the homogenised model of order $O(\varepsilon^3)$ in red (MSB-Stiff higher-order) and the finite element model in black. From left to right: (a) modal deformation in displacement $v(x)$, (b) difference between homogenised model and finite element model $v_{\text{Macro}}(x) - v_{\text{FE}}$, (c) modal deformation in rotation $\theta(x)$, and (d) difference between homogenised model and finite element model $\theta_{\text{Macro}}(x) - \theta_{\text{FE}}(x)$. The dots represent the nodes between the beams where the springs are located.</figcaption></figure>

### 4.3.3. Synthesis

In order to quantify the differences between the modes calculated according to a homogenised model and those calculated according to a finite element analysis, two indicators have been defined: (i) following (46), the modal deformations in displacement or in rotation are compared only at the nodes, such as the more $\overline{\text{MAC}}_{X=v \text{ or } \theta}$ is weak and
### 4.3.2. "Full" Macroscopic Modal Analysis

The full macroscopic modal analysis was carried out on a structure consisting of $N = 30$ periods. The two models used are that of the generic beam and that of the beam MSB-Stiff enriched with the higher orders; see Section 4.4. In each case, the four boundary conditions are those indicated in Section 4.2.

Figure 6 shows the fundamental modes determined by these two models with boundary conditions either fixed-fixed (Figure 6-top) or fixed-free (Figure 6-bottom). For comparison, the results obtained using finite elements and those obtained using "reduced" modal analysis are also shown.

Compared to the 10-period fixed-fixed structure, we can see that tripling the number of periods, as assumed, results in better scale separation and a gain in accuracy of the reduced modal analysis approach. Moreover, the concentration of the error for rotations on the periods at the edges is even more obvious.

For both types of boundary conditions, it can be seen that "full" modal analyses performed with four boundary conditions significantly improve the results. Furthermore, it appears that the higher-order model leads to a better approximation than the generic model. This is explained by the fact that the additional term in the higher-order model is inherited from effects specific to the contrast under study. This improves the description as it provides a model of an accuracy up to $O(\varepsilon^3)$. In contrast, the generic model is built up in such a way that it gives an accuracy up to $O(\varepsilon)$ regardless of the contrasts. It is always valid at the dominant order in the sense that the terms naturally vanish when they are not relevant. Nevertheless, even if the differential operator of the generic model is of fourth order, one cannot expect that it will reach the same level of accuracy as the higher-order model.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** First modal deformation of periodic structures ($N = 30$ periods, ($\beta = 10^{-2}$; $\delta = 1$)) with fixed-fixed (top) and fixed-free (bottom) boundary conditions. The homogenised models of order $O(\varepsilon)$ are shown in blue (MSB-Stiff) and green (MSB-Generic); the homogenised model of order $O(\varepsilon^3)$ in red (MSB-Stiff higher-order) and the finite element model in black. From left to right: (a) modal deformation in displacement $v(x)$, (b) difference between homogenised model and finite element model $v_{\text{Macro}}(x) - v_{\text{FE}}$, (c) modal deformation in rotation $\theta(x)$, and (d) difference between homogenised model and finite element model $\theta_{\text{Macro}}(x) - \theta_{\text{FE}}(x)$. The dots represent the nodes between the beams where the springs are located.</figcaption></figure>

### 4.3.3. Synthesis

In order to quantify the differences between the modes calculated according to a homogenised model and those calculated according to a finite element analysis, two indicators have been defined: (i) following (46), the modal deformations in displacement or in rotation are compared only at the nodes, such as the more $\overline{\text{MAC}}_{X=v \text{ or } \theta}$ is weak and
the more the compared modal shapes are close; (ii) the modal angular frequencies are compared according to (47), where $\gamma$ represents the local bending wavenumber; see (5).

$$ \widetilde{\mathrm{MAC}}_X = 1 - \mathrm{MAC}_X = 1 - \left( \frac{\underline{\underline{X}}^{\mathrm{homo}} \cdot \underline{\underline{X}}^{\mathrm{FE}}}{\|\underline{\underline{X}}^{\mathrm{homo}}\| \|\underline{\underline{X}}^{\mathrm{FE}}\|} \right)^2 \quad \underline{\underline{X}} \in \{\underline{v}, \theta\} \quad (46) $$

$$ \%\gamma = 100 \frac{|\gamma^{\mathrm{homo}} - \gamma^{\mathrm{FE}}|}{\gamma^{\mathrm{FE}}} \quad ; \quad \gamma^4 = \left(\frac{\omega}{\omega_b}\right)^2 \quad (47) $$

Figure 7 shows the three indicators $(\widetilde{\mathrm{MAC}}_v, \widetilde{\mathrm{MAC}}_\theta, \%\gamma)$ for the three first modes of the structure constituted by $N = 30$ periods. Systematically, it can be noted that (i) the “full” boundary conditions lead to more accurate descriptions, the best being obtained for the higher-order model, and (ii) the higher modes are described with less precision due to the smaller scale separation. It should be noted, however, that the dynamics of the system are relatively undisturbed by the “reduced” boundary conditions, as indicated by the good estimation of the eigenfrequencies ($\%\gamma \le 0.8\%$ for the third mode); see Figure 7-right.

<figure><img src="image_10.png" /><figcaption>Figure 7. Comparisons of the modes of the structures with $N = 30$ periods according to the indicators (from left to right) $(\widetilde{\mathrm{MAC}}_v, \widetilde{\mathrm{MAC}}_\theta, \%\gamma)$ for two types of boundary conditions: fixed-fixed (top), fixed-free (bottom). The homogenised models of order $O(\varepsilon)$ are shown in blue (MSB-Stiff) and green (MSB-Generic); the homogenised model of order $O(\varepsilon^3)$ in red (MSB-Stiff higher-order).</figcaption></figure>

### 4.4. Higher-Order Macroscopic Model

The aim is to enrich the MSB-Stiff model with higher order terms so that the order of the macroscopic differential operator corresponds to the order of the differential operator of the local problem (4) and allows all local conditions to be translated. In the enriched representation of discrete media, e.g., [39,41], and in higher-order continuous media, e.g., [25,26,31,42], it is necessary to construct generalised boundary conditions that are consistent with the order of the differential operator. On the other hand, for MSB-Stiff beams, the macroscopic model must be enriched to take into account all the real boundary conditions.

Returning to the asymptotic expansion of Section 3.2.1, solving it to second-order leads to the macroscopic kinematic relation (18b). Moreover, the transverse force equilibrium at order 2 gave the macroscopic balance at the dominant order. This equilibrium is rewritten in the form (48), which is systematically used in higher-order expressions.

$$ \frac{\rho L^2}{Ee_0^2} (\omega^{(1)})^2 \dot{v}^{(0)} = \beta^{(2)} \dot{v}^{(0)} - L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \partial_x^2 \dot{v}^{(0)} \quad (48) $$
the more the compared modal shapes are close; (ii) the modal angular frequencies are compared according to (47), where $\gamma$ represents the local bending wavenumber; see (5).

$$ \widetilde{\mathrm{MAC}}_X = 1 - \mathrm{MAC}_X = 1 - \left( \frac{\underline{\underline{X}}^{\mathrm{homo}} \cdot \underline{\underline{X}}^{\mathrm{FE}}}{\|\underline{\underline{X}}^{\mathrm{homo}}\| \|\underline{\underline{X}}^{\mathrm{FE}}\|} \right)^2 \quad \underline{\underline{X}} \in \{\underline{v}, \theta\} \quad (46) $$

$$ \%\gamma = 100 \frac{|\gamma^{\mathrm{homo}} - \gamma^{\mathrm{FE}}|}{\gamma^{\mathrm{FE}}} \quad ; \quad \gamma^4 = \left(\frac{\omega}{\omega_b}\right)^2 \quad (47) $$

Figure 7 shows the three indicators ($\widetilde{\mathrm{MAC}}_v, \widetilde{\mathrm{MAC}}_\theta, \%\gamma$) for the three first modes of the structure constituted by $N = 30$ periods. Systematically, it can be noted that (i) the “full” boundary conditions lead to more accurate descriptions, the best being obtained for the higher-order model, and (ii) the higher modes are described with less precision due to the smaller scale separation. It should be noted, however, that the dynamics of the system are relatively undisturbed by the “reduced” boundary conditions, as indicated by the good estimation of the eigenfrequencies ($\%\gamma \le 0.8\%$ for the third mode); see Figure 7-right.

<figure><img src="image_11.png" /><figcaption>Figure 7. Comparisons of the modes of the structures with $N = 30$ periods according to the indicators (from left to right) ($\widetilde{\mathrm{MAC}}_v, \widetilde{\mathrm{MAC}}_\theta, \%\gamma$) for two types of boundary conditions: fixed-fixed (top), fixed-free (bottom). The homogenised models of order $O(\varepsilon)$ are shown in blue (MSB-Stiff) and green (MSB-Generic); the homogenised model of order $O(\varepsilon^3)$ in red (MSB-Stiff higher-order).</figcaption></figure>

### 4.4. Higher-Order Macroscopic Model

The aim is to enrich the MSB-Stiff model with higher order terms so that the order of the macroscopic differential operator corresponds to the order of the differential operator of the local problem (4) and allows all local conditions to be translated. In the enriched representation of discrete media, e.g., [39,41], and in higher-order continuous media, e.g., [25,26,31,42], it is necessary to construct generalised boundary conditions that are consistent with the order of the differential operator. On the other hand, for MSB-Stiff beams, the macroscopic model must be enriched to take into account all the real boundary conditions.

Returning to the asymptotic expansion of Section 3.2.1, solving it to second-order leads to the macroscopic kinematic relation (18b). Moreover, the transverse force equilibrium at order 2 gave the macroscopic balance at the dominant order. This equilibrium is rewritten in the form (48), which is systematically used in higher-order expressions.

$$ \frac{\rho L^2}{Ee_0^2} (\omega^{(1)})^2 \dot{v}^{(0)} = \beta^{(2)} \dot{v}^{(0)} - L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \partial_x^2 \dot{v}^{(0)} \quad (48) $$
### 4.4.1. First-Order Macroscopic Model

The second-order moment (A4) gives a kinematic relation at the first order identical to that of the dominant order:

$$ \theta^{(1)} = \frac{1}{\delta^{(0)} + 1} \partial_x v^{(1)} \qquad (49) $$

The third-order shear equilibrium (A5), in which the relation (49) has been introduced, provides the first-order macroscopic equation:

$$ -L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \partial_x^2 v^{(1)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(1)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} v^{(0)} = 0 \qquad (50) $$

Denoting $\Theta^{(1)} = \theta^{(0)} + \epsilon\theta^{(1)}$ and $V^{(1)} = v^{(0)} + \epsilon v^{(1)}$ the first-order approximations of $\theta$ and $v$, respectively, the kinematic relation is written as (17) + $\epsilon$ (49):

$$ \Theta^{(1)} = \frac{1}{\delta + 1} \partial_x V^{(1)} $$

The asymptotic expansion of the angular frequency is $\omega = \epsilon\omega^{(1)} + \epsilon^2\omega^{(2)} + o(\epsilon^2)$. This involves that $\omega^2 = \epsilon^2(\omega^{(1)})^2 + 2\epsilon^3\omega^{(1)}\omega^{(2)} + o(\epsilon^3)$. Consequently, the first-order macroscopic equation is written as:

$$ m(\omega^2 - \omega_c^2) V^{(1)} + C_{eq} \partial_x^2 V^{(1)} = O(\epsilon^2) \qquad (51) $$

The Equation (51) is identical to (18a); only the precision of the expansion differs. Since the derivation order (two) of the operator is unchanged, it is not possible to apply all the boundary conditions. It is therefore necessary to expand asymptotically to the next order.

### 4.4.2. Second-Order Macroscopic Model

The second-order moment (A5) gives the second-order kinematic relation, rewritten using (48):

$$ \theta^{(2)} = \frac{1}{\delta^{(0)} + 1} \left( \partial_x v^{(2)} + \frac{\beta^{(2)}}{210} \left( 13 + \frac{1}{1 + \delta^{(0)}} \right) \partial_x v^{(0)} + \frac{L^2}{210} \frac{\delta^{(0)}}{\delta^{(0)} + 1} \left( 22 - \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right) \partial_x^3 v^{(0)} \right) \qquad (52) $$

The second-order approximations of $\theta$ and $v$ are denoted by

$$ \Theta^{(2)} = \theta^{(0)} + \epsilon\theta^{(1)} + \epsilon^2\theta^{(2)} \quad ; \quad V^{(2)} = v^{(0)} + \epsilon v^{(1)} + \epsilon^2 v^{(2)} $$

By summing (17) + $\epsilon$ (49) + $\epsilon^2$ (52), the kinematic relation enriched to the second order is:

$$ \Theta^{(2)} = \frac{C_{eq}}{C} \left( \partial_x V^{(2)} \left[ 1 + \frac{\beta}{210} \left( 13 + \frac{1}{1+\delta} \right) \right] + \frac{\ell^2}{210} \frac{\delta}{\delta+1} \left( 22 - \frac{\delta}{\delta+1} \right) \partial_x^3 V^{(2)} \right) + O(\epsilon^3) \qquad (53) $$

The fourth-order shear equilibrium (A6) provides the following second-order equilibrium equation, in which (48) has been used to swap the time and space derivatives:

$$ \beta^{(2)} v^{(2)} - \frac{1}{60} (\beta^{(2)})^2 v^{(0)} - \frac{\rho L^2}{E e_0^2} \left( (\omega^{(1)})^2 v^{(2)} + 2\omega^{(1)} \omega^{(2)} v^{(1)} + ((\omega^{(2)})^2 + 2\omega^{(1)} \omega^{(3)}) v^{(0)} \right) \\ - L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \left( \partial_x^2 v^{(2)} + \frac{\beta^{(2)}}{210} \left( 20 + \frac{1}{\delta^{(0)} + 1} \right) \partial_x^2 v^{(0)} \right) + \frac{L^4}{12} \left( 1 - \frac{2}{35} \left( \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right)^2 \left( \frac{23}{2} - \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right) \right) \partial_x^4 v^{(0)} = 0 \qquad (54) $$
### 4.4.1. First-Order Macroscopic Model

The second-order moment (A4) gives a kinematic relation at the first order identical to that of the dominant order:

$$ \theta^{(1)} = \frac{1}{\delta^{(0)} + 1} \partial_x v^{(1)} \qquad (49) $$

The third-order shear equilibrium (A5), in which the relation (49) has been introduced, provides the first-order macroscopic equation:

$$ -L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \partial_x^2 v^{(1)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(1)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} v^{(0)} = 0 \qquad (50) $$

Denoting $\Theta^{(1)} = \theta^{(0)} + \epsilon\theta^{(1)}$ and $V^{(1)} = v^{(0)} + \epsilon v^{(1)}$ the first-order approximations of $\theta$ and $v$, respectively, the kinematic relation is written as (17) + $\epsilon$ (49):

$$ \Theta^{(1)} = \frac{1}{\delta + 1} \partial_x V^{(1)} $$

The asymptotic expansion of the angular frequency is $\omega = \epsilon\omega^{(1)} + \epsilon^2\omega^{(2)} + o(\epsilon^2)$. This involves that $\omega^2 = \epsilon^2(\omega^{(1)})^2 + 2\epsilon^3\omega^{(1)}\omega^{(2)} + o(\epsilon^3)$. Consequently, the first-order macroscopic equation is written as:

$$ m(\omega^2 - \omega_c^2) V^{(1)} + C_{eq} \partial_x^2 V^{(1)} = O(\epsilon^2) \qquad (51) $$

The Equation (51) is identical to (18a); only the precision of the expansion differs. Since the derivation order (two) of the operator is unchanged, it is not possible to apply all the boundary conditions. It is therefore necessary to expand asymptotically to the next order.

### 4.4.2. Second-Order Macroscopic Model

The second-order moment (A5) gives the second-order kinematic relation, rewritten using (48):

$$ \theta^{(2)} = \frac{1}{\delta^{(0)} + 1} \left( \partial_x v^{(2)} + \frac{\beta^{(2)}}{210} \left( 13 + \frac{1}{1 + \delta^{(0)}} \right) \partial_x v^{(0)} + \frac{L^2}{210} \frac{\delta^{(0)}}{\delta^{(0)} + 1} \left( 22 - \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right) \partial_x^3 v^{(0)} \right) \qquad (52) $$

The second-order approximations of $\theta$ and $v$ are denoted by

$$ \Theta^{(2)} = \theta^{(0)} + \epsilon\theta^{(1)} + \epsilon^2\theta^{(2)} \quad ; \quad V^{(2)} = v^{(0)} + \epsilon v^{(1)} + \epsilon^2 v^{(2)} $$

By summing (17) + $\epsilon$ (49) + $\epsilon^2$ (52), the kinematic relation enriched to the second order is:

$$ \Theta^{(2)} = \frac{C_{eq}}{C} \left( \partial_x V^{(2)} \left[ 1 + \frac{\beta}{210} \left( 13 + \frac{1}{1+\delta} \right) \right] + \frac{\ell^2}{210} \frac{\delta}{\delta+1} \left( 22 - \frac{\delta}{\delta+1} \right) \partial_x^3 V^{(2)} \right) + O(\epsilon^3) \qquad (53) $$

The fourth-order shear equilibrium (A6) provides the following second-order equilibrium equation, in which (48) has been used to swap the time and space derivatives:

$$ \beta^{(2)} v^{(2)} - \frac{1}{60} (\beta^{(2)})^2 v^{(0)} - \frac{\rho L^2}{E e_0^2} \left( (\omega^{(1)})^2 v^{(2)} + 2\omega^{(1)} \omega^{(2)} v^{(1)} + ((\omega^{(2)})^2 + 2\omega^{(1)} \omega^{(3)}) v^{(0)} \right) - L^2 \frac{\delta^{(0)}}{\delta^{(0)} + 1} \left( \partial_x^2 v^{(2)} + \frac{\beta^{(2)}}{210} \left( 20 + \frac{1}{\delta^{(0)} + 1} \right) \partial_x^2 v^{(0)} \right) + \frac{L^4}{12} \left( 1 - \frac{2}{35} \left( \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right)^2 \left( \frac{23}{2} - \frac{\delta^{(0)}}{\delta^{(0)} + 1} \right) \right) \partial_x^4 v^{(0)} = 0 \qquad (54) $$
The sum of the three Equations (18a) + ε(50) + ε²(54) provides the second-order macroscopic behaviour equation in ε called MSB-Stiff higher-order:

$$m(\omega^2 - \omega_{c,eq}^2)V^{(2)}(x) + C^{(2)}\partial_x^2 V^{(2)}(x) - E^{(2)}I\ell\partial_x^4 V^{(2)}(x) = O(\epsilon^3) \quad (55)$$

with (C(2), E(2)) the elastic coefficients (C, E) enriched by taking into account the stiffnesses (β, δ).

$$C^{(2)} = C_{eq} \left[ 1 + \frac{\beta}{210} \left( 20 + \frac{1}{1+\delta} \right) \right] ; \quad E^{(2)} = E \left[ 1 - \frac{1}{35} \left( \frac{\delta}{1+\delta} \right)^2 \left( 23 - 2 \frac{\delta}{1+\delta} \right) \right] \quad (56)$$

The model (55) is enriched by two orders compared to MSB-Stiff (18a). This can be interpreted mechanically in terms of transverse forces and moments. This leads to the definition of beam model forces, generalised to the second order T(2)(x), M(2)(x), which satisfy the following system of equations:

$$\partial_x T^{(2)}(x) = \rho A \omega^2 V^{(2)}(x) - \frac{K_{eq}}{\ell} V^{(2)}(x) \quad (57a)$$

$$\partial_x M^{(2)}(x) + T^{(2)}(x) = -\frac{C^{(2)}}{\ell} \partial_x V^{(2)}(x) \quad (57b)$$

$$M^{(2)}(x) = -E^{(2)} I \partial_x V^{(2)}(x) \quad (57c)$$

$$\frac{C^{(2)}}{\ell} \partial_x V^{(2)}(x) = \frac{C_{eq} \left( 1 + \frac{\beta}{210} \left( 20 + \frac{1}{1+\delta} \right) \right)}{\ell \left( 1 + \frac{\beta}{210} \left( 13 + \frac{1}{1+\delta} \right) \right)} \left( \frac{C}{C_{eq}} \Theta^{(2)}(x) - \frac{\ell^2}{210} \frac{\delta}{\delta+1} \left( 22 - \frac{\delta}{\delta+1} \right) \partial_x^3 V^{(2)} \right) \quad (57d)$$

Indeed, it takes into account a bending moment M(2)(x), which raises the order of the differential operator to the fourth order. Note that the effective elastic modulus E(2) is modified by the presence of the rotation springs. The force balance (57a) involves a distributed force source whose elastic parameter K_{eq} takes into account both the compression spring support and the bending stiffness of the beam. Finally, the momentum balance (57b) is subjected to a source of distributed momentum associated with the gradient of the deflection. The elastic parameter C(2) involves both the compression and rotation spring support as well as the bending stiffness of the beam. The enriched kinematic relation shows that the deflection depends on the rotation Θ(2) and on the variation of curvature ∂_x^3 V(2).

### 4.4.3. Dispersion Curve Based on the Model MSB-Stiff Higher-Order

The dispersion equation of the model MSB-Stiff higher-order is obtained by inserting V(2) in the form V(2) = V*ε^ik_ωx into (55).

$$\left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta} \frac{C^{(2)}}{C_{eq}} (k_{\omega}\ell)^2 + \frac{E^{(2)}}{E} (k_{\omega}\ell)^4 \quad (58)$$

Figure 8 compares the dispersion relationships corresponding to the following models: (i) MSB-Stiff (34), (ii) MSB-Generic (37), (iii) MSB-Stiff higher-order (58), and (iv) Floquet-Bloch (39).

The MSB-Stiff higher-order model, like MSB-Generic, takes into account the bending moment and therefore the $(k_\omega \ell)^4$ term, which makes it possible to follow the variation in the dispersion branch more precisely than MSB-Stiff. For MSB-Stiff higher-order, the improvement over MSB-Generic lies in the introduction of correctors (C(2), E(2)), which are specific to the magnitudes of the contrasts $(\beta^{(2)}, \delta^{(0)})$. Improving the dispersion by integrating higher orders in the homogenised description is a well-known method; see, for example, [31] for composite elastic media, [43] for mass/spring networks, and [39] for 1D stratified media.
The sum of the three Equations (18a) + ε(50) + ε²(54) provides the second-order macroscopic behaviour equation in ε called MSB-Stiff higher-order:

$$m(\omega^2 - \omega_{c,eq}^2) V^{(2)}(x) + C^{(2)} \partial_x^2 V^{(2)}(x) - E^{(2)} I \ell \partial_x^4 V^{(2)}(x) = O(\epsilon^3) \quad (55)$$

with (C(2), E(2)) the elastic coefficients (C, E) enriched by taking into account the stiffnesses (β, δ).

$$C^{(2)} = C_{eq} \left[ 1 + \frac{\beta}{210} \left( 20 + \frac{1}{1+\delta} \right) \right] \quad ; \quad E^{(2)} = E \left[ 1 - \frac{1}{35} \left( \frac{\delta}{1+\delta} \right)^2 \left( 23 - 2 \frac{\delta}{1+\delta} \right) \right] \quad (56)$$

The model (55) is enriched by two orders compared to MSB-Stiff (18a). This can be interpreted mechanically in terms of transverse forces and moments. This leads to the definition of beam model forces, generalised to the second order T(2)(x), M(2)(x), which satisfy the following system of equations:

$$\partial_x T^{(2)}(x) = \rho A \omega^2 V^{(2)}(x) - \frac{K_{eq}}{\ell} V^{(2)}(x) \quad (57a)$$

$$\partial_x M^{(2)}(x) + T^{(2)}(x) = -\frac{C^{(2)}}{\ell} \partial_x V^{(2)}(x) \quad (57b)$$

$$M^{(2)}(x) = -E^{(2)} I \partial_x V^{(2)}(x) \quad (57c)$$

$$\frac{C^{(2)}}{\ell} \partial_x V^{(2)}(x) = \frac{C_{eq}}{\ell} \frac{1 + \frac{\beta}{210} \left(20 + \frac{1}{1+\delta}\right)}{1 + \frac{\beta}{210} \left(13 + \frac{1}{1+\delta}\right)} \left( \frac{C}{C_{eq}} \Theta^{(2)}(x) - \frac{\ell^2}{210} \frac{\delta}{\delta+1} \left( 22 - \frac{\delta}{\delta+1} \right) \partial_x^3 V^{(2)} \right) \quad (57d)$$

Indeed, it takes into account a bending moment M(2)(x), which raises the order of the differential operator to the fourth order. Note that the effective elastic modulus E(2) is modified by the presence of the rotation springs. The force balance (57a) involves a distributed force source whose elastic parameter K_{eq} takes into account both the compression spring support and the bending stiffness of the beam. Finally, the momentum balance (57b) is subjected to a source of distributed momentum associated with the gradient of the deflection. The elastic parameter C(2) involves both the compression and rotation spring support as well as the bending stiffness of the beam. The enriched kinematic relation shows that the deflection depends on the rotation Θ(2) and on the variation of curvature ∂_x^3 V(2).

### 4.4.3. Dispersion Curve Based on the Model MSB-Stiff Higher-Order

The dispersion equation of the model MSB-Stiff higher-order is obtained by inserting V(2) in the form V(2) = V* e ik_ω x into (55).

$$\left(\frac{\omega}{\omega_b}\right)^2 = \left(\frac{\omega_{c,eq}}{\omega_b}\right)^2 + \frac{12\delta}{1+\delta} \frac{C^{(2)}}{C_{eq}} (k_{\omega} \ell)^2 + \frac{E^{(2)}}{E} (k_{\omega} \ell)^4 \quad (58)$$

Figure 8 compares the dispersion relationships corresponding to the following models: (i) MSB-Stiff (34), (ii) MSB-Generic (37), (iii) MSB-Stiff higher-order (58), and (iv) Floquet-Bloch (39).

The MSB-Stiff higher-order model, like MSB-Generic, takes into account the bending moment and therefore the $(k_\omega \ell)^4$ term, which makes it possible to follow the variation in the dispersion branch more precisely than MSB-Stiff. For MSB-Stiff higher-order, the improvement over MSB-Generic lies in the introduction of correctors (C(2), E(2)), which are specific to the magnitudes of the contrasts $(\beta^{(2)}, \delta^{(0)})$. Improving the dispersion by integrating higher orders in the homogenised description is a well-known method; see, for example, [31] for composite elastic media, [43] for mass/spring networks, and [39] for 1D stratified media.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Dispersion relationships (large wavelength) for the system of contrasts ($p = 2, q = 0$) with $\epsilon = 0.1$ according to four models: MSB-Stiff (solid line), MSB-Generic (dotted line), MSB-Stiff higher-order (dashed line), and Floquet-Bloch (dots).</figcaption></figure>

**5. Conclusions**

The discrete medium homogenization method (HPDM) applied to multi-supported beams yields three distinct models (29)–(31), which correspond to stiff, medium, and soft rotation supports, respectively. Independently of the rotation springs, the stiffness of the compression springs gives rise to the emergence of a cut-off angular frequency below which the macroscopic dynamics is not propagative. The different models are synthesised into a single generic model that fits the Floquet-Bloch dispersion curve up to wavelengths of the order of three periods. Furthermore, the homogenization also provides an interpretation in terms of macroscopic forces/moments.

A mixed “discrete/continuous” approach to the boundary conditions is proposed, which allows the boundary conditions (natural or essential) actually applied at the local scale to be expressed in terms of Robin-type boundary conditions on macroscopic variables. The modal analysis performed with these boundary conditions and the generic model gives results in good agreement with a full finite element calculation, with great economy of numerical resources. Note that this condition can be used to define interface conditions that allow different models to be connected.

Nevertheless, for sufficiently rigid supports whose homogenised model at the dominant order is of differential order two, it is shown that the results provided by the generic model and the “discrete/continuous” type boundary conditions can be further improved by the dedicated homogenised model of higher order.

This work will be expanded in several directions. The findings of this study can be extended to Timoshenko beams. The developments will be of a similar nature but will require the introduction of another dimensionless parameter concerning the relative effects of shear versus bending in the beams. This additional parameter will influence the macroscopic model without modifying the general framework of the models or the treatment of the boundary conditions.

Furthermore, it is straightforward to integrate a linear dissipative contribution to the purely conservative supports in order to address, for instance, visco-elastic polymer fixation systems or to deal with support foundations on soils expressed by impedance conditions. In order to consider such supports, it is sufficient to consider springs with complex frequency-dependent values. In this case, the dissipation introduced is transferred to the macroscopic model, whose apparent dissipation and its frequency dependence will depend on the beam/support system.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Dispersion relationships (large wavelength) for the system of contrasts ($p = 2, q = 0$) with $\varepsilon = 0.1$ according to four models: MSB-Stiff (solid line), MSB-Generic (dotted line), MSB-Stiff higher-order (dashed line), and Floquet-Bloch (dots).</figcaption></figure>

**5. Conclusions**

The discrete medium homogenization method (HPDM) applied to multi-supported beams yields three distinct models (29)–(31), which correspond to stiff, medium, and soft rotation supports, respectively. Independently of the rotation springs, the stiffness of the compression springs gives rise to the emergence of a cut-off angular frequency below which the macroscopic dynamics is not propagative. The different models are synthesised into a single generic model that fits the Floquet-Bloch dispersion curve up to wavelengths of the order of three periods. Furthermore, the homogenization also provides an interpretation in terms of macroscopic forces/moments.

A mixed “discrete/continuous” approach to the boundary conditions is proposed, which allows the boundary conditions (natural or essential) actually applied at the local scale to be expressed in terms of Robin-type boundary conditions on macroscopic variables. The modal analysis performed with these boundary conditions and the generic model gives results in good agreement with a full finite element calculation, with great economy of numerical resources. Note that this condition can be used to define interface conditions that allow different models to be connected.

Nevertheless, for sufficiently rigid supports whose homogenised model at the dominant order is of differential order two, it is shown that the results provided by the generic model and the “discrete/continuous” type boundary conditions can be further improved by the dedicated homogenised model of higher order.

This work will be expanded in several directions. The findings of this study can be extended to Timoshenko beams. The developments will be of a similar nature but will require the introduction of another dimensionless parameter concerning the relative effects of shear versus bending in the beams. This additional parameter will influence the macroscopic model without modifying the general framework of the models or the treatment of the boundary conditions.

Furthermore, it is straightforward to integrate a linear dissipative contribution to the purely conservative supports in order to address, for instance, visco-elastic polymer fixation systems or to deal with support foundations on soils expressed by impedance conditions. In order to consider such supports, it is sufficient to consider springs with complex frequency-dependent values. In this case, the dissipation introduced is transferred to the macroscopic model, whose apparent dissipation and its frequency dependence will depend on the beam/support system.
It is also worth noting the possibility of dealing with high frequencies, which results in the modulation of local beam modes following the approach implemented in [44,45]. The transposition of the mixed "discrete/continuous" boundary conditions to these high-frequency dynamic regimes makes it possible to envisage a better consideration of real boundary conditions.

The same method could also be used to analyse the instability of multi-supported beams by considering the static equilibrium of bended beams in their deformed state under the action of the normal force instead of the dynamic equilibrium expressed on the undeformed state as assumed in this study.

Finally, these diverse perspectives also extend to the analysis of the mechanical behaviour of multi-supported plates or two-dimensional networks of multi-supported beams at their nodes.

**Author Contributions:** Conceptualization, A.R. and C.B.; methodology, A.R. and C.B.; formal analysis, A.R. and C.B.; writing—original draft preparation, A.R. and C.B.; writing—review and editing, A.R. and C.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Data are contained within the article.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

FE Finite Element

HPDM Homogenization of Periodic and Discrete Media

MSB Multi-Supported Beam

## Appendix A. Expansion of Equations of Equilibrium by Order of ε

In this appendix are presented the full expansions of the transverse and momentum equilibria (16) to successive first orders of ε:

### Order -1:

$$ \begin{cases} T^{(-1)} = 0 \\ M^{(-1)} = \delta^{(-1)} L \theta^{(-1)} = 0 \end{cases} \qquad (A1) $$

### Order 0:

$$ \begin{cases} T^{(0)} = 0 \\ M^{(0)} = \delta^{(0)} L \theta^{(0)} + (\delta^{(0)} + 1) L \theta^{(-1)} = 0 \end{cases} \qquad (A2) $$

### Order 1:

$$ \begin{cases} T^{(1)} = L^2 \partial_x \theta^{(-1)} = 0 \\ M^{(1)} = \delta^{(-1)} L \theta^{(1)} + (\delta^{(0)} + 1) L \theta^{(0)} + \delta^{(1)} L \theta^{(-1)} - L \partial_x v^{(0)} = 0 \end{cases} \qquad (A3) $$

### Order 2:

$$ \begin{cases} T^{(2)} = L^2 \partial_x \theta^{(0)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} - L^2 \partial_x^2 v^{(0)} = 0 \\ M^{(2)} = \delta^{(-1)} L \theta^{(2)} + (\delta^{(0)} + 1) L \theta^{(1)} + \delta^{(1)} L \theta^{(0)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210 E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(-1)} - L \partial_x v^{(1)} + \frac{L^3}{6} \partial_x^2 \theta^{(-1)} \end{cases} \qquad (A4) $$

### Order 3:
It is also worth noting the possibility of dealing with high frequencies, which results in the modulation of local beam modes following the approach implemented in [44,45]. The transposition of the mixed “discrete/continuous” boundary conditions to these high-frequency dynamic regimes makes it possible to envisage a better consideration of real boundary conditions.

The same method could also be used to analyse the instability of multi-supported beams by considering the static equilibrium of bended beams in their deformed state under the action of the normal force instead of the dynamic equilibrium expressed on the undeformed state as assumed in this study.

Finally, these diverse perspectives also extend to the analysis of the mechanical behaviour of multi-supported plates or two-dimensional networks of multi-supported beams at their nodes.

**Author Contributions:** Conceptualization, A.R. and C.B.; methodology, A.R. and C.B.; formal analysis, A.R. and C.B.; writing—original draft preparation, A.R. and C.B.; writing—review and editing, A.R. and C.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Data are contained within the article.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

FE Finite Element

HPDM Homogenization of Periodic and Discrete Media

MSB Multi-Supported Beam

## Appendix A. Expansion of Equations of Equilibrium by Order of ε

In this appendix are presented the full expansions of the transverse and momentum equilibria (16) to successive first orders of ε:

### Order -1:

$$ \begin{cases} T^{(-1)} = 0 \\ M^{(-1)} = \delta^{(-1)} L \theta^{(-1)} = 0 \end{cases} \qquad (A1) $$

### Order 0:

$$ \begin{cases} T^{(0)} = 0 \\ M^{(0)} = \delta^{(0)} L \theta^{(0)} + (\delta^{(0)} + 1) L \theta^{(-1)} = 0 \end{cases} \qquad (A2) $$

### Order 1:

$$ \begin{cases} T^{(1)} = L^2 \partial_x \theta^{(-1)} = 0 \\ M^{(1)} = \delta^{(-1)} L \theta^{(1)} + (\delta^{(0)} + 1) L \theta^{(0)} + \delta^{(1)} L \theta^{(-1)} - L \partial_x v^{(0)} = 0 \end{cases} \qquad (A3) $$

### Order 2:

$$ \begin{cases} T^{(2)} = L^2 \partial_x \theta^{(0)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(0)} - L^2 \partial_x^2 v^{(0)} = 0 \\ M^{(2)} = \delta^{(-1)} L \theta^{(2)} + (\delta^{(0)} + 1) L \theta^{(1)} + \delta^{(1)} L \theta^{(0)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210 E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(-1)} - L \partial_x v^{(1)} + \frac{L^3}{6} \partial_x^2 \theta^{(-1)} \end{cases} \qquad (A4) $$

### Order 3:
$$
\left\{
\begin{aligned}
\mathbb{T}^{(3)} &= L^2 \partial_x \theta^{(1)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(1)} - L^2 \partial_x^2 v^{(1)} \\
&\quad + \left[ \beta^{(3)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] v^{(0)} + \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L^2 \partial_x \theta^{(-1)} + \frac{L^4}{6} \partial_x^3 \theta^{(-1)} = 0 \\
\mathbb{M}^{(3)} &= \delta^{(-1)} L \theta^{(3)} + (\delta^{(0)} + 1) L \theta^{(2)} + \delta^{(1)} L \theta^{(1)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(0)} - L \partial_x v^{(2)} + \frac{L^3}{6} \partial_x^2 \theta^{(0)} \\
&\quad + \left[ \delta^{(3)} - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] L \theta^{(-1)} - \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L \partial_x v^{(0)} - \frac{L^3}{6} \partial_x^3 v^{(0)}
\end{aligned}
\right.
\quad (A5)
$$

Order 4:

$$
\left\{
\begin{aligned}
\mathbb{T}^{(4)} &= L^2 \partial_x \theta^{(2)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(2)} - L^2 \partial_x^2 v^{(2)} + \left[ \beta^{(3)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] v^{(1)} \\
&\quad + \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L^2 \partial_x \theta^{(0)} + \frac{L^4}{6} \partial_x^3 \theta^{(0)} + \left[ \beta^{(4)} - \frac{\rho L^2}{E e_0^2} (\omega^{(2)})^2 - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(3)} - \left( \frac{\rho L^2}{E e_0^2} \right)^2 \frac{(\omega^{(1)})^4}{60} \right] v^{(0)} \\
&\quad - \frac{9\rho L^2}{70E e_0^2} (\omega^{(1)})^2 L^2 \partial_x^2 v^{(0)} - \frac{L^4}{12} \partial_x^4 v^{(0)} + \frac{13\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} L^2 \partial_x^2 \theta^{(-1)} = 0 \\
\mathbb{M}^{(4)} &= \delta^{(-1)} L \theta^{(4)} + (\delta^{(0)} + 1) L \theta^{(3)} + \delta^{(1)} L \theta^{(2)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(1)} - L \partial_x v^{(3)} + \frac{L^3}{6} \partial_x^2 \theta^{(1)} \\
&\quad + \left[ \delta^{(3)} - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] L \theta^{(0)} - \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L \partial_x v^{(1)} - \frac{L^3}{6} \partial_x^3 v^{(1)} \\
&\quad + \left[ \delta^{(4)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(2)})^2 - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(3)} - \frac{13}{970200} \left( \frac{\rho L^2}{E e_0^2} \right)^2 (\omega^{(1)})^4 \right] L \theta^{(-1)} \\
&\quad + \frac{\rho L^2}{140E e_0^2} (\omega^{(1)})^2 L^3 \partial_x^2 \theta^{(-1)} - \frac{13\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} L \partial_x v^{(0)} + \frac{L^5}{72} \partial_x^4 \theta^{(-1)} = 0
\end{aligned}
\right.
\quad (A6)
$$

References

1. Colmenares, D.; Andersson, A.; Karoumi, R. Closed-form solution for mode superposition analysis of continuous beams on flexible supports under moving harmonic loads. *J. Sound Vib.* **2022**, *520*, 116587. [CrossRef]

2. Li, J.; Chen, B.; Mao, H. Exact closed-form solution for vibration characteristics of multi-span beams on an elastic foundation subjected to axial force. *Structures* **2024**, *60*, 105884. [CrossRef]

3. Mathews, P.M. Vibrations of a beam on elastic foundation. *ZAMM J. Appl. Math. Mech. Z. Angew. Math. Mech.* **1958**, *38*, 105–115. [CrossRef]

4. Achenbach, J.D.; Sun, C.T. Moving load on a flexibly supported Timoshenko beam. *Int. J. Solids Struct.* **1965**, *1*, 353–370. [CrossRef]

5. Chonan, S. Moving Harmonic Load on an Elastically Supported Timoshenko Beam. *ZAMM J. Appl. Math. Mech. Z. Angew. Math. Mech.* **1978**, *58*, 9–15. [CrossRef]

6. Metrikine, A.V.; Verichev, S.N. Instability of vibrations of a moving two-mass oscillator on a flexibly supported Timoshenko beam. *Arch. Appl. Mech.* **2001**, *71*, 613–624. [CrossRef]

7. Rodrigues, C.; Simões, F.M.; Pinto da Costa, A.; Froio, D.; Rizzi, E. Finite element dynamic analysis of beams on nonlinear elastic foundations under a moving oscillator. *Eur. J. Mech.-A/Solids* **2018**, *68*, 9–24. [CrossRef]

8. Morfidis, K. Vibration of Timoshenko beams on three-parameter elastic foundation. *Comput. Struct.* **2010**, *88*, 294–308. [CrossRef]

9. Ding, L.; Wu, L.; Zhu, H.P. Propagation and localization of wave in multi-span Timoshenko beams on elastic foundations under moving harmonic loads. *Int. J. Numer. Anal. Methods Geomech.* **2017**, *41*, 1687–1710. [CrossRef]

10. Mead, D.J. A new method of analyzing wave propagation in periodic structures; Applications to periodic timoshenko beams and stiffened plates. *J. Sound Vib.* **1986**, *104*, 9–27. [CrossRef]
$$
\left\{
\begin{aligned}
T^{(3)} &= L^2 \partial_x \theta^{(1)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(1)} - L^2 \partial_x^2 v^{(1)} \\
&\quad + \left[ \beta^{(3)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] v^{(0)} + \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L^2 \partial_x \theta^{(-1)} + \frac{L^4}{6} \partial_x^3 \theta^{(-1)} = 0 \\
M^{(3)} &= \delta^{(-1)} L \theta^{(3)} + (\delta^{(0)} + 1) L \theta^{(2)} + \delta^{(1)} L \theta^{(1)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(0)} - L \partial_x v^{(2)} + \frac{L^3}{6} \partial_x^2 \theta^{(0)} \\
&\quad + \left[ \delta^{(3)} - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] L \theta^{(-1)} - \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L \partial_x v^{(0)} - \frac{L^3}{6} \partial_x^3 v^{(0)}
\end{aligned}
\right.
\quad (A5)
$$

Order 4:

$$
\left\{
\begin{aligned}
T^{(4)} &= L^2 \partial_x \theta^{(2)} + \left[ \beta^{(2)} - \frac{\rho L^2}{E e_0^2} (\omega^{(1)})^2 \right] v^{(2)} - L^2 \partial_x^2 v^{(2)} + \left[ \beta^{(3)} - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] v^{(1)} \\
&\quad + \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L^2 \partial_x \theta^{(0)} + \frac{L^4}{6} \partial_x^3 \theta^{(0)} + \left[ \beta^{(4)} - \frac{\rho L^2}{E e_0^2} (\omega^{(2)})^2 - \frac{\rho L^2}{E e_0^2} 2\omega^{(1)} \omega^{(3)} - \left( \frac{\rho L^2}{E e_0^2} \right)^2 \frac{(\omega^{(1)})^4}{60} \right] v^{(0)} \\
&\quad - \frac{9\rho L^2}{70E e_0^2} (\omega^{(1)})^2 L^2 \partial_x^2 v^{(0)} - \frac{L^4}{12} \partial_x^4 v^{(0)} + \frac{13\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} L^2 \partial_x^2 \theta^{(-1)} = 0 \\
M^{(4)} &= \delta^{(-1)} L \theta^{(4)} + (\delta^{(0)} + 1) L \theta^{(3)} + \delta^{(1)} L \theta^{(2)} + \left[ \delta^{(2)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(1)})^2 \right] L \theta^{(1)} - L \partial_x v^{(3)} + \frac{L^3}{6} \partial_x^2 \theta^{(1)} \\
&\quad + \left[ \delta^{(3)} - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} \right] L \theta^{(0)} - \frac{13\rho L^2}{210E e_0^2} (\omega^{(1)})^2 L \partial_x v^{(1)} - \frac{L^3}{6} \partial_x^3 v^{(1)} \\
&\quad + \left[ \delta^{(4)} - \frac{\rho L^2}{210E e_0^2} (\omega^{(2)})^2 - \frac{\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(3)} - \frac{13}{970200} \left( \frac{\rho L^2}{E e_0^2} \right)^2 (\omega^{(1)})^4 \right] L \theta^{(-1)} \\
&\quad + \frac{\rho L^2}{140E e_0^2} (\omega^{(1)})^2 L^3 \partial_x^2 \theta^{(-1)} - \frac{13\rho L^2}{210E e_0^2} 2\omega^{(1)} \omega^{(2)} L \partial_x v^{(0)} + \frac{L^5}{72} \partial_x^4 \theta^{(-1)} = 0
\end{aligned}
\right.
\quad (A6)
$$

References

1. Colmenares, D.; Andersson, A.; Karoumi, R. Closed-form solution for mode superposition analysis of continuous beams on flexible supports under moving harmonic loads. *J. Sound Vib.* **2022**, 520, 116587. [CrossRef]

2. Li, J.; Chen, B.; Mao, H. Exact closed-form solution for vibration characteristics of multi-span beams on an elastic foundation subjected to axial force. *Structures* **2024**, 60, 105884. [CrossRef]

3. Mathews, P.M. Vibrations of a beam on elastic foundation. *ZAMM J. Appl. Math. Mech. Z. Angew. Math. Mech.* **1958**, 38, 105–115. [CrossRef]

4. Achenbach, J.D.; Sun, C.T. Moving load on a flexibly supported Timoshenko beam. *Int. J. Solids Struct.* **1965**, 1, 353–370. [CrossRef]

5. Chonan, S. Moving Harmonic Load on an Elastically Supported Timoshenko Beam. *ZAMM J. Appl. Math. Mech. Z. Angew. Math. Mech.* **1978**, 58, 9–15. [CrossRef]

6. Metrikine, A.V.; Verichev, S.N. Instability of vibrations of a moving two-mass oscillator on a flexibly supported Timoshenko beam. *Arch. Appl. Mech.* **2001**, 71, 613–624. [CrossRef]

7. Rodrigues, C.; Simões, F.M.; Pinto da Costa, A.; Froio, D.; Rizzi, E. Finite element dynamic analysis of beams on nonlinear elastic foundations under a moving oscillator. *Eur. J. Mech.-A/Solids* **2018**, 68, 9–24. [CrossRef]

8. Morfidis, K. Vibration of Timoshenko beams on three-parameter elastic foundation. *Comput. Struct.* **2010**, 88, 294–308. [CrossRef]

9. Ding, L.; Wu, L.; Zhu, H.P. Propagation and localization of wave in multi-span Timoshenko beams on elastic foundations under moving harmonic loads. *Int. J. Numer. Anal. Methods Geomech.* **2017**, 41, 1687–1710. [CrossRef]

10. Mead, D.J. A new method of analyzing wave propagation in periodic structures; Applications to periodic timoshenko beams and stiffened plates. *J. Sound Vib.* **1986**, 104, 9–27. [CrossRef]
11. Mead, D.J.; Yaman, Y. The harmonic response of uniform beams on multiple linear supports: A flexural wave analysis. *J. Sound Vib.* **1990**, *141*, 465–484. [CrossRef]

12. Lin, Y.K.; McDaniel, T.J. Dynamics of Beam-Type Periodic Structures. *J. Eng. Ind. 1969*, *91*, 1133–1141.

13. Lin, H.P.; Chang, S.C. Free vibration analysis of multi-span beams with intermediate flexible constraints. *J. Sound Vib.* **2005**, *281*, 155–169. [CrossRef]

14. Manevitch, L.I.; Oshmyan, V.G. An Asymptotic Study of the Linear Vibrations of a Stretched Beam with Concentrated Masses and Discrete Elastic Supports. *J. Sound Vib.* **1999**, *223*, 679–691. [CrossRef]

15. Sanchez-Palencia, E. *Non-Homogeneous Media and Vibration Theory*; Springer: Berlin/Heidelberg, Geramny, 1980. [CrossRef]

16. Auriault, J.L.; Boutin, C.; Geindreau, C. *Homogenization of Coupled Phenomena in Heterogenous Media*; Wiley-ISTE: London, UK; Hoboken, NJ, USA, 2009. [CrossRef]

17. Allaire, G. Homogenization and Two-Scale Convergence. *SIAM J. Math. Anal.* **1992**, *23*, 1482–1518. [CrossRef]

18. Caillerie, D.; Trompette, P.; Verna, P. Homogenization of periodic trusses. In Proceedings of the Congress IASS Madrid, 1989; pp. 7139–7180.

19. Moreau, G.; Caillerie, D. Continuum modeling of lattice structures in large displacement applications to buckling analysis. *Comput. Struct.* **1998**, *68*, 181–189. [CrossRef]

20. Tollenaere, H.; Caillerie, D. Continuous modeling of lattice structures by homogenization. *Adv. Eng. Softw.* **1998**, *29*, 699–705. [CrossRef]

21. Boutin, C.; Hans, S.; Chesnais, C. Generalized Beams and Continua. Dynamics of Reticulated Structures. In *Advances in Mechanics and Mathematics*; Springer: New York, NY, USA, 2010; Volume 21, pp. 131–141. [CrossRef]

22. Hans, S.; Boutin, C. Dynamics of discrete framed structures: A unified homogenized description. *J. Mech. Mater. Struct* **2008**, *3*, 1709–1739. [CrossRef]

23. Kolpakov, A.G.; Andrianov, I.V. Asymptotic decomposition in the problem of joined elastic beams. *ZAMM Z. Angew. Math. Mech.* **2014**, *94*, 818–836. [CrossRef]

24. Gambin, B.; Kröner, E. Higher-Order Terms in the Homogenized Stress-Strain Relation of Periodic Elastic Media. *Phys. Status Solidi* (b) **1989**, *151*, 513–519. [CrossRef]

25. Germain, P. The Method of Virtual Power in Continuum Mechanics. Part 2: Microstructure. *SIAM J. Appl. Math.* **1973**, *25*, 556–575. [CrossRef]

26. Mindlin, R.D. Second gradient of strain and surface-tension in linear elasticity. *Int. J. Solids Struct.* **1965**, *1*, 417–438. [CrossRef]

27. Eringen, A.C. Linear theory of nonlocal elasticity and dispersion of plane waves. *Int. J. Eng. Sci.* **1972**, *10*, 425–435. [CrossRef]

28. Boutin, C. Microstructural effects in elastic composites. *Int. J. Solids Struct.* **1996**, *33*, 1023–1051. [CrossRef]

29. Smyshlyaev, V.P.; Cherednichenko, K.D. On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media. *J. Mech. Phys. Solids* **2000**, *48*, 1325–1357. [CrossRef]

30. Cherednichenko, K.; Smyshlyaev, V.P. On full two-scale expansion of the solutions of nonlinear periodic rapidly oscillating problems and higher-order homogenised variational problems. *Arch. Ration. Mech. Anal.* **2004**, *174*, 385–442. [CrossRef]

31. Boutin, C.; Auriault, J. Rayleigh scattering in elastic composite materials. *Int. J. Eng. Sci.* **1993**, *31*, 1669–1689. [CrossRef]

32. Fish, J.; Chen, W. Higher-Order Homogenization of Initial/Boundary-Value Problem. *J. Eng. Mech.* **2001**, *127*, 1223–1230. [CrossRef]

33. Andrianov, I.V.; Bolshakov, V.I.; Danishevs'kyy, V.V.; Weichert, D. Higher order asymptotic homogenization and wave propagation in periodic composite materials. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2008**, *464*, 1181–1201. [CrossRef]

34. Dumontet, H. Study of a boundary layer problem in elastic composite materials. *ESAIM Math. Model. Numer. Anal.* **1986**, *20*, 265–286. [CrossRef]

35. Moskow, S.; Vogelius, M. First-order corrections to the homogenised eigenvalues of a periodic composite medium. A convergence proof. *Proc. R. Soc. Edinb. Sect. A Math.* **1997**, *127*, 1263–1299. [CrossRef]

36. Allaire, G.; Amar, M. Boundary layer tails in periodic homogenization. *ESAIM Control Optim. Calc. Var.* **1999**, *4*, 209–243. [CrossRef]

37. Cioranescu, D.; Donato, P. *An Introduction to Homogenization*; Oxford University Press: Oxford, UK, 1999. [CrossRef]

38. Armstrong, S.; Kuusi, T.; Mourrat, J.C.; Prange, C. Quantitative Analysis of Boundary Layers in Periodic Homogenization. *Arch. Ration. Mech. Anal.* **2017**, *226*, 695–741. [CrossRef]

39. Cornaggia, R.; Guzina, B.B. Second-order homogenization of boundary and transmission conditions for one-dimensional waves in periodic media. *Int. J. Solids Struct.* **2020**, *188–189*, 88–102. [CrossRef]

40. Fergoug, M.; Parret-Fréaud, A.; Feld, N.; Marchand, B.; Forest, S. A general boundary layer corrector for the asymptotic homogenization of elastic linear composite structures. *Compos. Struct.* **2022**, *285*, 115091. [CrossRef]

41. Kaplunov, J.D.; Pichugin, A.V. *On Rational Boundary Conditions for Higher-Order Long-Wave Models*; IUTAM Bookseries; Springer: Dordrecht, The Netherlands, 2009, Volume 10, pp. 81–90. [CrossRef] [PubMed]

42. Andrianov, I.V.; Bulanova, N.S.; Sedin, V.L. Vibrations of ribbed plates on elastic bases. *Int. Appl. Mech.* **1999**, *35*, 64–68. [CrossRef]

43. Pichugin, A.V.; Askes, H.; Tyas, A. Asymptotic equivalence of homogenisation procedures and fine-tuning of continuum theories. *J. Sound Vib.* **2008**, *313*, 858–874. [CrossRef]
11. Mead, D.J.; Yaman, Y. The harmonic response of uniform beams on multiple linear supports: A flexural wave analysis. *J. Sound Vib.* **1990**, *141*, 465–484. [CrossRef]

12. Lin, Y.K.; McDaniel, T.J. Dynamics of Beam-Type Periodic Structures. *J. Eng. Ind. 1969*, *91*, 1133–1141.

13. Lin, H.P.; Chang, S.C. Free vibration analysis of multi-span beams with intermediate flexible constraints. *J. Sound Vib.* **2005**, *281*, 155–169. [CrossRef]

14. Manevitch, L.I.; Oshmyan, V.G. An Asymptotic Study of the Linear Vibrations of a Stretched Beam with Concentrated Masses and Discrete Elastic Supports. *J. Sound Vib.* **1999**, *223*, 679–691. [CrossRef]

15. Sanchez-Palencia, E. *Non-Homogeneous Media and Vibration Theory*; Springer: Berlin/Heidelberg, Geramny, 1980. [CrossRef]

16. Auriault, J.L.; Boutin, C.; Geindreau, C. *Homogenization of Coupled Phenomena in Heterogenous Media*; Wiley-ISTE: London, UK; Hoboken, NJ, USA, 2009. [CrossRef]

17. Allaire, G. Homogenization and Two-Scale Convergence. *SIAM J. Math. Anal.* **1992**, *23*, 1482–1518. [CrossRef]

18. Caillerie, D.; Trompette, P.; Verna, P. Homogenization of periodic trusses. In Proceedings of the Congress IASS Madrid, 1989; pp. 7139–7180.

19. Moreau, G.; Caillerie, D. Continuum modeling of lattice structures in large displacement applications to buckling analysis. *Comput. Struct.* **1998**, *68*, 181–189. [CrossRef]

20. Tollenaere, H.; Caillerie, D. Continuous modeling of lattice structures by homogenization. *Adv. Eng. Softw.* **1998**, *29*, 699–705. [CrossRef]

21. Boutin, C.; Hans, S.; Chesnais, C. Generalized Beams and Continua. Dynamics of Reticulated Structures. In *Advances in Mechanics and Mathematics*; Springer: New York, NY, USA, 2010; Volume 21, pp. 131–141. [CrossRef]

22. Hans, S.; Boutin, C. Dynamics of discrete framed structures: A unified homogenized description. *J. Mech. Mater. Struct* **2008**, *3*, 1709–1739. [CrossRef]

23. Kolpakov, A.G.; Andrianov, I.V. Asymptotic decomposition in the problem of joined elastic beams. *ZAMM Z. Angew. Math. Mech.* **2014**, *94*, 818–836. [CrossRef]

24. Gambin, B.; Kröner, E. Higher-Order Terms in the Homogenized Stress-Strain Relation of Periodic Elastic Media. *Phys. Status Solidi* (b) **1989**, *151*, 513–519. [CrossRef]

25. Germain, P. The Method of Virtual Power in Continuum Mechanics. Part 2: Microstructure. *SIAM J. Appl. Math.* **1973**, *25*, 556–575. [CrossRef]

26. Mindlin, R.D. Second gradient of strain and surface-tension in linear elasticity. *Int. J. Solids Struct.* **1965**, *1*, 417–438. [CrossRef]

27. Eringen, A.C. Linear theory of nonlocal elasticity and dispersion of plane waves. *Int. J. Eng. Sci.* **1972**, *10*, 425–435. [CrossRef]

28. Boutin, C. Microstructural effects in elastic composites. *Int. J. Solids Struct.* **1996**, *33*, 1023–1051. [CrossRef]

29. Smyshlyaev, V.P.; Cherednichenko, K.D. On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media. *J. Mech. Phys. Solids* **2000**, *48*, 1325–1357. [CrossRef]

30. Cherednichenko, K.; Smyshlyaev, V.P. On full two-scale expansion of the solutions of nonlinear periodic rapidly oscillating problems and higher-order homogenised variational problems. *Arch. Ration. Mech. Anal.* **2004**, *174*, 385–442. [CrossRef]

31. Boutin, C.; Auriault, J. Rayleigh scattering in elastic composite materials. *Int. J. Eng. Sci.* **1993**, *31*, 1669–1689. [CrossRef]

32. Fish, J.; Chen, W. Higher-Order Homogenization of Initial/Boundary-Value Problem. *J. Eng. Mech.* **2001**, *127*, 1223–1230. [CrossRef]

33. Andrianov, I.V.; Bolshakov, V.I.; Danishevs'kyy, V.V.; Weichert, D. Higher order asymptotic homogenization and wave propagation in periodic composite materials. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2008**, *464*, 1181–1201. [CrossRef]

34. Dumontet, H. Study of a boundary layer problem in elastic composite materials. *ESAIM Math. Model. Numer. Anal.* **1986**, *20*, 265–286. [CrossRef]

35. Moskow, S.; Vogelius, M. First-order corrections to the homogenised eigenvalues of a periodic composite medium. A convergence proof. *Proc. R. Soc. Edinb. Sect. A Math.* **1997**, *127*, 1263–1299. [CrossRef]

36. Allaire, G.; Amar, M. Boundary layer tails in periodic homogenization. *ESAIM Control Optim. Calc. Var.* **1999**, *4*, 209–243. [CrossRef]

37. Cioranescu, D.; Donato, P. *An Introduction to Homogenization*; Oxford University Press: Oxford, UK, 1999. [CrossRef]

38. Armstrong, S.; Kuusi, T.; Mourrat, J.C.; Prange, C. Quantitative Analysis of Boundary Layers in Periodic Homogenization. *Arch. Ration. Mech. Anal.* **2017**, *226*, 695–741. [CrossRef]

39. Cornaggia, R.; Guzina, B.B. Second-order homogenization of boundary and transmission conditions for one-dimensional waves in periodic media. *Int. J. Solids Struct.* **2020**, *188–189*, 88–102. [CrossRef]

40. Fergoug, M.; Parret-Fréaud, A.; Feld, N.; Marchand, B.; Forest, S. A general boundary layer corrector for the asymptotic homogenization of elastic linear composite structures. *Compos. Struct.* **2022**, *285*, 115091. [CrossRef]

41. Kaplunov, J.D.; Pichugin, A.V. *On Rational Boundary Conditions for Higher-Order Long-Wave Models*; IUTAM Bookseries; Springer: Dordrecht, The Netherlands, 2009, Volume 10, pp. 81–90. [CrossRef] [PubMed]

42. Andrianov, I.V.; Bulanova, N.S.; Sedin, V.L. Vibrations of ribbed plates on elastic bases. *Int. Appl. Mech.* **1999**, *35*, 64–68. [CrossRef]

43. Pichugin, A.V.; Askes, H.; Tyas, A. Asymptotic equivalence of homogenisation procedures and fine-tuning of continuum theories. *J. Sound Vib.* **2008**, *313*, 858–874. [CrossRef]
44. Craster, R.V.; Kaplunov, J.; Pichugin, A.V. High-frequency homogenization for periodic media. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2010**, *466*, 2341–2362. [CrossRef]

45. Rallu, A.; Hans, S.; Boutin, C. Asymptotic analysis of high-frequency modulation in periodic systems. Analytical study of discrete and continuous structures. *J. Mech. Phys. Solids* **2018**, *117*, 123–156. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
44. Craster, R.V.; Kaplunov, J.; Pichugin, A.V. High-frequency homogenization for periodic media. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2010**, *466*, 2341–2362. [CrossRef]

45. Rallu, A.; Hans, S.; Boutin, C. Asymptotic analysis of high-frequency modulation in periodic systems. Analytical study of discrete and continuous structures. *J. Mech. Phys. Solids* **2018**, *117*, 123–156. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.