Computational study of a homogenized nonlinear generalization of Timoshenko beam proposed by Turco et al.

1.   Introduction

Because of the increasing demands in engineering applications ^[20,47], the need for advanced engineering materials is increasing, too. Mechanical metamaterials are becoming relevant in this regard in that they can achieve high strength-to-density ratio, vibration control, and energy absorption ^[42] properties. Mechanical metamaterials owe their performances to their relevant micro-structure–at a given sub-length-scale–rather than to the inherent characteristics of the manufacturing materials ^[45]. Mechanical metamaterials are currently being studied for potential application to civil engineering, aeronautical engineering, mechanical engineering–especially soft robotics–and biomedical engineering ^[24,37,50]. Mechanical metamaterials have been proved as viable for protecting infrastructures against strong explosions and impacts and, more generally, damage ^[19,21,52], thanks to their excellent energy-absorption capabilities ^[7,26]. Still concerning civil engineering, they have also been explored to mitigate the consequences of ordinary impacts that may occur in accidents, e.g., in highways ^[50]. Thanks to their high strength-to-density ratios, they could also be used as resistant and lightweight materials enabling the redesign of some traditional civil infrastructure systems ^[6]. Metamaterials can be designed so as to be capable of controlling passive vibrations, which confers to them the potential to improve the seismic resilience of buildings ^[51] and bridges ^[40]. While not strictly mechanical, metamaterials with enhanced hygro-thermo-mechanical properties may also be employed for building's façade design, to achieve low energy consumption and improve indoor comfort ^[29].

Beam theory is of great importance for the development of efficient mathematical models to be used for obtaining design predictions, whose empirical test would be not viable economically or technically. In most cases, when dealing with mechanical metamaterials, one deals with large deformations and displacements, for which nonlinear models must be employed ^[33]. The small-scale slender structural elements that make up lattice mechanical metamaterials, simply called sometimes microbeams, are also susceptible to large deformations ^{[32,35,36,48,49]}. Therefore, the use of nonlinear models for describing these elements is necessary for the prediction of the macroscopic mechanical phenomena involving metamaterials ^{[1,16,30,31,41]}. Moreover, when the slenderness is insufficiently high, it is necessary to consider shear deformation effects, which are neglected in the classical Euler-Bernoulli theory ^[11,27]. The Timoshenko model can consider these effects, although it is originally a model for small deformations ^[2,28]. In that sense, classical beam theories of continuum media cannot be applied for the description of metamaterial microstructures, so it is necessary to enrich such theories ^{[8,9,12,17,22,38]}. Timoshenko's theory, however, has stood out for its multiple applications in various engineering fields where a more accurate analysis of beam deformations is required. Therefore, many discretizations have been formulated to be able to apply it in multiple contexts, being a valuable inspiration for the creation of generalized models ^[43]. In the work by Giorgio et al. ^[23], a 2D continuous model was compared with a 3D microstructured model, considering in both the geometry of pantograph networks. The simulations evidenced not only that the 2D model is computationally less expensive, but also that satisfactory and similar results were produced between both models. This was achieved despite visually ignoring the details of each beam element. Despite this, working with mesoscale discrete models, element by element, can easily become complex, especially for periodically repeated structural elements such as metamaterials. For that reason, it may be more convenient to work with homogenized continuum models of discrete models, known as multi-scale modeling ^[3].

Multiscale modeling, however, is a great challenge due to multiple reasons. First of all, complex interactions at small scales produce exotic behavior at the macroscopic scale of the metamaterial that does not necessarily follow the physical principles that are fulfilled at the small scale. However, it is this exotic behavior that most often is the target of metamaterial design, so accurate modeling of this coupling between scales is necessary. As a consequence, this coupling is very difficult to model accurately using general continuum models that consider the system as a whole, i.e., direct continuum formulations ^[10]. Discrete models, along with computational methods, can help in understanding this coupling, at the usual cost of demanding computational needs. Asymptotic homogenization techniques can exploit the understanding brought by discrete models to get continuum models that properly take into account the coupling between scales ^[34] in metamaterials and, hence, to get an understanding of the metamaterial at the macroscopic level ^[46]. Asymptotic homogenization techniques are, however, mostly heuristic and, even when rigorous convergence results are available, there is no guarantee that the obtained homogenized continuum model adequately captures the desired behaviors of the original discrete model, since asymptotic homogenization intrinsically implies an information loss ^[5]. Indeed, continuum models offer a simplified version of the discrete model by averaging the microscopic scale behaviors.

The nonlinear discrete elastic model formulated by Turco et al. ^[43] is introduced as a generalization of the Timoshenko model. This model is indeed able to consider not only bending deformation effects, but also axial and shear deformation effects, which could not necessarily be neglected in mechanical metamaterials. It is discrete in nature and can analyze large deformations in beams, hence proving suitable for the design of mechanical metamaterials. For example, Eremeyev and Turco ^[18] employed it in the numerical modeling of pantograph-shaped metamaterial beams. Furthermore, Barchiesi et al. ^[4] studied a duoskelion (Greek word meaning 'two legs') beam making use of this discrete model. Although it was possible to satisfactorily predict in several benchmarks or more application-oriented cases the equilibrium shapes of beams at large deformations with the discrete model introduced in ^[43], its continuum counterpart, obtained by heuristic asymptotic homogenization, must still be validated against the discrete model. This may allow for the identification of possible issues and the assessment of the reliability of the model for design. Thus, the question arises: How are large deformations in beams manifested using the continuum homogenized nonlinear Timoshenko model for mechanical metamaterials? Aimed at addressing this question, the present numerical study, performed using the commercial finite element software COMSOL Multiphysics, takes into account several loading cases and elucidates the sensitivity of the homogenized continuum with respect to axial, bending, and shear stiffness parameters, as well as to load imperfections, in terms of mechanical behavior, including buckling onset and post-critical behavior. The main novelty of this paper lies in the numerical computation of the continuum homogenized model derived as the asymptotic limit, when the cell size tends to zero, of a nonlinear generalization of the Timoshenko beam proposed by Turco et al in ^[43], with the objective of validating its predictions with respect to those of the discrete model. The plan of the work is the following. In Section 2, we recall the employed discrete model and its homogenization. In Section 3, we report the results of some tests, taken from ^[43], performed by means of the homogenized model: cantilever beam subjected to varying transverse end-load, cantilever beam subjected to varying axial end-load with varying imperfection, clamped-hinged semicircular arch subjected to varying concentrated pressure force at its midpoint, hinged-hinged semicircular arch subjected to concentrated pressure force at its midpoint, and cantilever beam subjected to varying axial-transverse end-load leading to a curled configuration. Section 4 reports the conclusions of the present study and future outlooks.

2.   Modeling

This section focuses on describing the mathematical formulation and homogenization of the model proposed by Turco et al. ^[43]. As mentioned in the introduction, it corresponds to a generalization of the Timoshenko-Ehrenfest beam; see the work by Elishakoff ^[13].

2.1.   Discrete model

For the discrete model, an approach in the spirit of Hencky ^[25] is considered. Indeed, having in mind the augmented kinematic of the Timoshenko beam, a discrete beam is introduced made up of several straight elements that may not nencessarily form a rectilinear line. These elements are connected at nodes, as shown in Figure 1. Each node has a position $P_j$ and $p_j$, where $j = 1, 2, ..., N$, in the reference and current configuration, respectively. In the reference configuration, each straight element is associated with a unit vector $\mathbf{D}_{1, j} = \frac{P_{j+1}-P_j}{\|P_j+1-P_j\|}$. A triad of unit vectors orthogonal to each other, namely, {$\mathbf{D}_{1, j}, \mathbf{D}_{2, j}, \mathbf{D}_{3, j}$}, is also considered, with $\mathbf{D}_{3, j} = \mathbf{D}_{1, j}\times\mathbf{D}_{2, j}$. This orthogonal triad transforms in the current configuration into {$\mathbf{d}_{1, j} = \mathbf{Q}_{j}\mathbf{D}_{1, j}, \mathbf{d}_{2, j} = \mathbf{Q}_{j}\mathbf{D}_{2, j}, \mathbf{d}_{3, j} = \mathbf{Q}_{j}\mathbf{D}_{3, j}$}, where $\mathbf{Q}_{j}$ is a rotation tensor. Referring to Figure 1, we define the deformation measure $\Delta\mathbf{w}_{j+1}$ through the following vector subtraction:

Axial deformation $\Delta\mathbf{w}_{a, j+1}$ and shear deformation $\Delta\mathbf{w}_{c, j+1}$ measures are defined for each straight element as, respectively:

and

where this last definition makes use of the deformation $\Delta\mathbf{w}_{j+1}$. We now seek to define the tensor $\Delta\mathbf{Q}_{j+1}$, which represents a measure of bending deformation of the discrete beam. It is defined as the relative rotation between the two straight elements concurring at node $j+1$. Restricting ourselves to in-plane motions, denoting with $\phi_{j}$ the rotation angle (positive when rotation occurs in the counterclockwise direction) associated to the tensor $\mathbf{Q}_{j}$, we are led to the following elementary matrix-form expression for the tensor $\mathbf{Q}_{j}$ in the othornormal basis $(\mathrm{\mathbf{e}}_{1}, \mathrm{\mathbf{e}}_{2})$ in Figure 1:

Having now expressed the tensor $\mathbf{Q}_{j}$ in terms of the rotation angle $\phi_{j}$, we define the bending measure $\Delta\mathbf{Q}_{j+1}$ of the beam as the difference between the relative rotations of the current configuration and the reference configuration, such that:

that in direct form reads as:

After having introduced the kinematic quantities of the discrete model describing the current configuration of the system, we can introduce the deformation energy of the discrete beam, which is defined as the sum over the whole system of the following elementary contributions, namely, the axial/extensional contribution:

the bending contribution:

and the shear contribution:

The first and the third are associated to a single straight element, while the second to two adjacent straight elements. The quantities $\Delta\mathbf{w}_{a, j+1}$ and $\Delta\mathbf{w}_{c, j+1}$ are vectorial strain measures, therefore, the energies are defined making use of their Euclidean norm. The quantity $\Delta\mathbf{Q}_{j+1}$ is a tensorial strain measure, its norm appears in the bending contribution, and its square is computed as the sum of the squares of each component:

which can be expressed in direct form as:

2.2.   Homogenization toward the continuum model

We now focus on the derivation of the continuum model describing the discrete one illustrated in the previous subsection. To this end, we first declare the independent Lagrangian parameters for the continuum model and their relationship with those of the discrete model. Thereafter, we will take the number of straight elements to infinity, keeping fixed the total length of the considered system. Hence, the length $||{P}_{j+1}-{P}_{j}|| = \ell$ of the straight elements will be let to zero.

Let us now define the arclength abscissa $s\in[0, L]$ on the curve $\chi_0:[0, L]\rightarrow \mathbb{E}^{2}$ with length $L$ that gathers all the points ${P}_{j}$ in the continuum model. The independent Lagrangian parameters of the continuum model are declared to be the functions:

The function $\chi$ positions the one-dimensional continuum beam in the two-dimensional space $\mathbb{E}^{2}$, and hence describes–from a microscopic standpoint–the motion of the points ${P}_{j}$. On the other hand, the function $\phi$ gives the rotation angle of the beam cross-section, and, hence–from a microscopic standpoint–describes the rotation of all unit vectors $\mathbf{D}_{1, j}$. Mathematically, the following relationships between the continuum independent Lagrangian parameters and the discrete ones are considered as:

It is worth mentioning that, henceforth, anytime differentiation will be performed, indicated by a prime, it will be done with respect to $s$. For convenience, the expression $\mathbf{Q}_{j}\mathbf{D}_{1, j} = \mathbf{d}_{1, j}$ will be denoted as $\mathbf{e}(\phi_j)$. The Lagrangian parameters $\chi$ and $\phi$ are expressed through Taylor expansions performed with respect to $s$ truncated at the first order as:

Using elementary trigonometric identities, the quantity $\|\Delta\mathbf{Q}_{j+1}\|^{2} = \mathrm{tr}(\Delta\mathbf{Q}_{j+1}^{T}\Delta\mathbf{Q}_{j+1})$ can be written as:

Using the Taylor expansion in Eq (2.15), the quantity $\cos(\phi(s_{j+1})-\phi(s_{j}))$ can be approximated at the second order as:

that can be replaced along with Eq (2.14) in the Eqs (2.7)–(2.9). With an abuse of notation, denoting with $\phi_j$ the angle formed by $\mathbf{d}_{1, j}$ with a fixed line, any line parallel to $\mathbf{e}_1$, as an instance, we obtain the deformation energy of the discrete system in terms of the Lagrangian parameters of the continuum computed at $s_j$:

where $\kappa_0$ is the initial curvature of the beam, namely that of the curve that is parametrized over the arclength abscissa $s$, which can be obtained as $\kappa_0 = \chi_0^{\prime\prime}\cdot\chi_0^\prime$. We now specify the scaling law for micro-scale stiffnesses, by defining the macro-scale stiffnesses $\underline{a}$, $\underline{b}$ and $\underline{c}$ which relate with the axial, bending and shear stiffness, respectively, of the discrete model as follows:

The scaling law above, where the macro-scale stiffnesses $\underline{a}$, $\underline{b}$, and $\underline{c}$ are independent of the cell size $\ell$, means that the micro-scale stiffnesses scale with $\ell$ in the same way, thus yielding in the limit of $\ell$ approaching zero a system where none of the stiffnesses is negligible with respect to the others. To get the continuum model, the limit $\ell\rightarrow0$ is considered, neglecting the higher order infinitesimals $o(\ell)$. Indeed, this is the error generated in the homogenization process when passing from the discrete model to the continuum. In conclusion, the deformation energy of the homogenized continuum model reads as:

It is noted that, the homogenized deformation energy above characterizes the continuum as one of micropolar type ^[14,15,44] and that, when $\nabla \mathrm{\mathbf{u}}\approx\mathbf{0}$ and $\phi\approx0$, it coincides with that of a classical Timoshenko beam. The deformation above energy is needed to compute equilibrium configurations predicted by the model. The weak form equation that is discretized and then solved numerically within the framework of a Galerkin approach is obtained by equating to zero the first variation of the potential energy, which is defined as the difference between the deformation energy and the work done by external generalized forces. In formulas, we have that the potential energy reads as

and, assuming that only external boundary forces and couples can be possibly acting onto the system, it can be re-written as

where $\mathrm{\mathbf{F}}_{s}$ and $M_{s}$, with $s\in\{0, L\}$, are boundary forces and couples ($M_{s}$ is positive when it induces a counterclockwise rotation), respectively, and $\mathrm{\mathbf{u}}(\cdot)$ is the displacement function defined as

where $\mathrm{\mathbf{e}}_{1}$ is the first unit vector of the orthonormal basis $(\mathrm{\mathbf{e}}_{1}, \mathrm{\mathbf{e}}_{2})$ in Figure 1. Replacing the displacement in the potential energy formula, taking into account the expression of the deformation energy, we get

where $\left(\cdot\right)_{1}$ and $\left(\cdot\right)_{2}$ denote the first and second components of the argument in the parentheses in the basis $(\mathrm{\mathbf{e}}_{1}, \mathrm{\mathbf{e}}_{2})$, respectively, and $u_{1}$ and $u_{2}$ are the first and second components of the displacement vector $\mathrm{\mathbf{u}}$ in the basis $(\mathrm{\mathbf{e}}_{1}, \mathrm{\mathbf{e}}_{2})$, respectively. After simplification, the expression above reads as

The Galerkin method is used to approximate the fields $u_{1}$, $u_{2}$ and $\phi$ appearing in the weak form. Specifically, each field has been approximated by considering its projection over the finite dimensional space of Lagrangian quadratic polynominal built over the spatial mesh. A simple equispaced mesh has been considered and the weak form has been inserted in the COMSOL software using the Weak Form PDE interface. As mentioned above, the weak form equation governing the system is obtained by equating to zero the first variation $\delta U$ of the potential energy, namely

where $\widetilde{U}$ is the integrand of the deformation energy and the two components of $\frac{\partial\widetilde{U}}{\partial\mathrm{\mathbf{u}}^{\prime}}$ are given according to the formula

and the formula

while the quantities $\frac{\partial\widetilde{U}}{\partial\phi}$ and $\frac{\partial\widetilde{U}}{\partial\phi^{\prime}}$ are given by the simple expressions $\frac{\partial\widetilde{U}}{\partial\phi} = c\left(\frac{\chi^{\prime}\cdot\mathrm{\mathbf{e}_{\perp}(\phi)}}{\Vert\chi^{\prime}\Vert}\right)$, where the simbol $\perp$ means that the vector $\mathbf{e}(\phi)$ has been rotated by an angle of $\pi/2$ in the anticlockwise direction, and $\frac{\partial\widetilde{U}}{\partial\phi^{\prime}} = 2b\left[\phi^{\prime}-\left(\chi_{0}^{\prime\prime}\right)_{1}\left(\chi_{0}^{\prime}\right)_{1}-\left(\chi_{0}^{\prime\prime}\right)_{2}\left(\chi_{0}^{\prime}\right)_{2}\right]$, respectively. Making use of a standard variational deduction procedure, in the Appendix A we provide the governing continuum equations for the studied beam as well as possible consistent boundary conditions.

3.   Numerical results

The main results obtained from the computational study of the beam model introduced in the previous section are reported in this section. The governing Euler-Lagrange equations have been numerically solved in weak form by means of the commercial finite element software COMSOL Multiphysics. Both straight and curved beams have been considered. The main acting force will be denoted with the symbol $F$, while the (relatively small) load imperfections will be denoted with the symbol $\varepsilon$. It is worth mentioning that, except for the last case that was carried out in two steps, the figures showing the evolution of the deformed shape with the loading step will make use of the blue color to indicate the unloaded/stress-free configuration, of the orange color to indicate the final deformed configuration, and of intermediate colors for correspondingly intermediate configurations. The bending stiffness $b$ is assumed to take the value $0.021$ N throughout the whole paper.

3.1.   Buckling of a cantilever beam

The problem of analyzing the pre- and post-critical buckling of a cantilever beam in large displacements and deformations has been addressed. The problem, including the analyzed loading conditions, is illustrated graphically in Figure 2. The beam has a length L equal to $1$ m and a dead increasing axial load F is applied at its free end. In addition, a small non-monotonously variable dead transverse load p(F) is applied, which, as the loading step increases, reaches a maximum value equal to $\varepsilon$F$_{f}$, where F$_{f}$ is the maximum applied axial force, which is reached at the final loading step. Figure 3 reports graphically the behavior of the imperfection load. The ratios among the stiffnesses of the beam have been chosen so as to fulfill $aL^2/b = 10000$ and $cL^2/b = 10000$, which implies a predominance of the bending behavior. It is to be noted that the previous ratios do have a physical dimension. Indeed, while they were considered in the paper by Turco et al. ^[43]–and this is why we are considering them here–to have nondimensional comparable quantities, due to the scaling laws and the appearance of the integration operator, they are no more nondimensional.

Simulation results for this test are reported in Figure 4 and Figure 5. Figure 4 reports the evolution of the deformed shape of the beam as the loading step increases, with the parameter $\varepsilon$ equal to $5\times10^{-4}$, while Figure 5 reports the force-displacement diagram for several values of the parameter $\varepsilon$. More specifically, up to scaling, it reports the force $F$ against the transverse displacement $u_2$ of the free end of the cantilever beam. It is seen that, as the imperfection load decreases, it affects in an increasingly less significant way the force-displacement diagram, while still keeping its beneficial effect in terms of observing the occurrence of buckling and, hence, compute the post-buckling behavior. Due to the higher axial stiffness coefficient–compared to the bending stiffness coefficient–the beam initially resists by mainly deforming only axially, until the imperfection, due to, e.g., load eccentricity, makes the transverse displacement sufficiently relevant for the axial load to further increase dramatically the non-axial deformation. The higher the ratio $a/b$ and the parameter $\varepsilon$, the sooner the buckling-like behavior onset is observed. Proper buckling behavior is observed when $\varepsilon$ becomes extremely small and it can be seen that the quantity $FL^2 / b$ in Figure 5, at buckling onset, which occurs in the small deformation linear regime, converges to the value $\pi^2 / 2$, which can be analytically related with the Euler elastic buckling load.

3.2.   Cantilever beam subjected to a transverse force at its free end

The objective of considering this case is to have a comparison with one of the cases reported in Turco et al. ^[43] in discrete form, where a cantilever beam subjected to a transverse end load is considered as shown in Figure 6. The beam has a length L equal to $1$ m and a dead increasing transverse load F is applied at its free end. Several stiffness ratios were considered in the analysis to assess the sensitivity of the beam model behavior upon the constitutive parameters: at first, the ratios $aL^2/b = 100$ and $cL^2/b = 10000$ were considered, then $aL^2/b = 10000$ and $cL^2/b = 10000$, and, finally, $aL^2/b = 10000$ and $cL^2/b = 10$.

Simulation results for this test are reported in Figures 7 and 8. Figure 7 reports the evolution of the deformed shape of the beam as the loading step increases for the considered stiffness ratios, while Figure 8 reports the force-displacement diagram for the considered stiffness ratios. More specifically, up to scaling, it reports the force $F$ against the transverse displacement $u_2$ of the free end of the cantilever beam. For the first scenario in Figure 7a, due to the high value of the axial and shear stiffnesses with respect to the bending stiffness, the beam tends to deform mostly by bending, initially. Subsequently and after having reached a certain deflection, the beam is positioned almost vertically and deforms axially until the end, which is why we observe a progressive stiffening (see Figure 8a). Figure 7b reports the same behavior, except for stiffening that is more pronounced due to the significantly higher value considered for the axial stifness with respect to the bending one. In the third scenario in Figure 7c, although again a similar behavior occurs, the curvature of the beam is mostly concentrated at its clamped end towards the final part of the test, the great majority of the beam being straight and, hence, not subjected to bending. When comparing the results of the first scenario with those obtained by Turco et al. ^[43]–second and third scenario were not analyzed therein–similar behaviors can be observed.

3.3.   Buckling of semicircular arches subjected to a vertical force at their midpoint

The cases analyzed in this sub-section concern semicircular arches in large deformation regime and seek to compare some simulations reported by Turco et al. ^[43] with the model in its discrete form. The studied problems, including the analyzed loading conditions, are illustrated graphically in Figure 9. The beam has a length L equal to $\pi$ m and a dead increasing vertical load F is applied at its midpoint. In addition, for the third analyzed problem Figure 9c, a small constant dead horizontal load $\varepsilon$ is applied, where $\varepsilon$ is equal to $0.01F_f$, being $F_f$ the final vertical force applied in the third scenario. The considered semicircular arches all have the same radius $r = 1$ m and, therefore, a uniform curvature $\kappa_0 = 1/r$. The three studied scenarios consider stiffness ratios such as $aL^2/b = 1000$ and $cL^2/b = 1000$, hence being bending predominant.

The simulation results are shown in Figures 10 and 11. For the first scenario in Figure 10a, it is observed that, when the right end of the arc reaches a certain horizontal displacement, said end tarts moving in the opposite direction. This is because the arch tends to maintain its total length due to its relatively high axial stiffness with respect to the bending stiffness. Different deformation modes can be found for the same arch, altering the boundary conditions and adding small forces as imperfections. The results are very similar to those reported in Turco et al. ^[43].

3.4.   Curled equilibrium configuration of a cantilever beam subjected to a transverse force

The simulation reported in this sub-section deals with the large deformation of a cantilever beam subjected to two sets of external forces that are applied in a sequence. The objective of this two-part loading procedure is to prove that there exist two different solution for the cantilever beam problem with an applied transverse end load. Such a loading procedure indeed allows for computing a curled solution when a purely transverse, i.e., vertical, force is applied. Clearly, when a purely transverse force is applied at the free end of a cantilever beam, one also has the more classical solution showing only deflection, see Figure 7a. The simulation presented in this sub-section aims at comparing with the results obtained by Turco et al. (2020) using the discrete form of the model, as well as experimental equilibrium configurations reported in Baroudi et al. (2019). The two sets of applied forces are visualized in Figure 12. In the first part of the loading procedure, an increasing axial compression force $\mu$ is applied along with a small constant transverse force $\varepsilon$ pointing upward, that constitutes an imperfection load. In the second part of the loading procedure, a decreasing axial tensile force $\mu$ is applied along with an increasing transverse vertical force $F$ pointing downward, until the force $\mu$ vanishes and one gets the desired solution. The stiffness ratios employed for this simulation are $aL^2/b = 100$ and $cL^2/b = 10000$, namely those leading to the results in Figure 7a. The same final vertical force utilized in the second test case is considered, too.

The simulation results of the test case analyzed in this sub-section are shown in Figure 13 and Figure 14. Results for the first part of the loading procedure are reported in blue, while those for the second one are reported in magenta. In the first loading part, the action of the constant imperfection force $\varepsilon = 0.003\mu_{f}$, being 0.003 an arbitrary value, allows sudden buckling along with the action of the force $\mu$. In the second part of the loading procedure the forces that act on the beam are linearly varied so as to being led at the final step to an acting vertical force $F$ only.

3.5.   Mesh convergence

In order to demonstrate that the number of mesh elements is sufficient for the simulations to give accurate results, the energy has been plotted against the load parameter, using different mesh elements, considering the test case of Figure 9b. Figure 15 reports such results and it is observed that the number of elements used for all simulations gives a solution that is adequately accurate with respect to that obtained by making use of exponentially finer meshes.

4.   Conclusions

The present investigation carried out the analysis of a homogenized continuum model that has been developed to model the large deformation mechanical regime of one-dimensional systems of particles having enriched kinematics, which is inspired by the Timoshenko model, considering not only bending deformations, but also axial and shear deformations. With the study of four test cases, we have tried to elucidate the features of the continuum and discrete models, as well as to assess to what extent the discrete model can be considered a good approximation of the continuum, and viceversa. At first, it shall be remarked that, although the results shown were mostly calculated with 1000 loading steps, calculations with up to 100 loading steps can be considered sufficient for obtaining accurate force-displacement diagrams, especially if there are no sudden changes in the curves, and can be obtained in a few seconds. By using the same stiffness ratios as those utilized in Turco et al. ^[43], the deformation results are quite similar, hence proving that the homogenized continuous model developed by Turco et al. ^[43] is suitable for metamaterials design and synthesis, being sufficiently descriptive of its discrete counterpart, hence being capable to model shear-type elements. At the same time, the continuum is more efficient than its discrete counterpart for performing computational calculations when the number of straight elements is very large. For small systems, the continuum is still comparable with the naturally discrete approach in terms of computational cost, since even with as few elements as about sixteen, the accuracy of the finite element solution of the continuum seems satisfactory. All considered, the continuum model analyzed in this paper may yield viable moving forward with the development of future research towards the evolution of two-dimensional or three-dimensional models for the analysis of complex mechanical metamaterials, which need to consider a greater number of elements in their structure.

We conclude this contribution by pointing out the fact that, although the homogenized model has given satisfactory results, it is still necessary to judiciously and reasonably establish the rigidities of the model and the applied loads. Very low values of the stiffnesses or incorrect relationships between them can cause the solution to diverge, as well as very high values of the loads or insufficient loading conditions. On the other hand, it may be beneficial for the model to study or incorporate more sophisticated calculation methods than Newton's method, such as the arc length method of ^[39], which can consider bifurcation behaviors and above limit forces. Furthermore, it is also advisable to carry out studies on the possibility of considering dynamic effects when, for example, the mass of the beam cannot be neglected.

Authors contribution

J.M.T.E. conceived the presented idea and performed the computations. E.B. conceived the presented idea, developed the theory, performed the computations, and supervised the research. All authors discussed the results and contributed to the final manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research reported in the present contribution was carried out as part of the project "Metamaterials design and synthesis with applications to infrastructure engineering" funded by the MUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022–grant 20228CPHN5.

Conflict of interest

Emilio Barchiesi is an editorial boardmember for [Networks and Heterogeneous Media] and was not involved inthe editorial review or the decision to publish this article. The authors declare there is no conflict of interest.

Appendix A

Making use of a standard variational deduction procedure, we here provide the governing continuum equations for the studied beam as well as possible consistent boundary conditions. Integrating by parts the weak form Eq (2.32), we get by standard arguments of calculus of variations the governing Euler-Lagrange equations and associated boundary conditions. The Euler-Lagrange equations read as

where

and the expression of $\frac{\mathrm{d}}{\mathrm{d}s}\frac{\partial\widetilde{U}}{\partial\mathrm{\mathbf{u}}^{\prime}}$ reads component-wise as

and

Finally, the boundary conditions read as