Modular modeling approach for FDM printed structures and piezo disks for metamaterial design

Research in metamaterials has recently gained interest in the field of noise and vibration control. The ability of creating band gap zones with minimum added mass is the main feature behind its success. In addition, the use of 3D printed parts, particularly the Fused Deposition Modeling (FDM), offers a practical solution for manufacturing parts with intricate shapes that are challenging for standard manufacturing processes, which is the case of many structural metamaterials. The combination of this concept with smart materials can further improve performance, providing the means to overcome typical issues. From a design perspective, the problem with coupling rises, as both the mechanical and electrical responses relies on the load circuit and on the mechanical properties of the smart elements. Therefore, the modeling of such structures is a rather complex task, for it involves multiphysical simulations of systems with complex geometries and typically a high number of degrees of freedom. Hence, this paper presents a direct approach for modeling and simulating smart metamaterials using a state-space formulation, which allows the modular coupling of electromechanical resonators manufactured by FDM. The numerical results are compared to experimental data obtained with unit cells prototypes embedded with piezoelectric elements and connected with a tunable shunt electric circuit. The good agreement between test and simulated data validates the design procedure.


INTRODUCTION
Research in metamaterials is gaining an ever-increasing interest in the field of noise and vibration control (Essink & Inman 2016;Miranda Júnior, Ferreira, & dos Santos 2017;Yang, Lee, & Kim 2016).The ability of creating frequency band gap zones, in which transfer of mechanical energy is limited, without penalizing other structural parameters, is the main purpose of this kind of solution.Firstly achieved in a mechanical system by Liu et al. (2000), the band gap mechanism stands on the Bragg scattering effect.Later, subwavelength band gap mechanisms were proposed based on local resonances and wave lensing (Cummer;Christensen;Alù, 2016;Huang;Shen;Jing, 2016).In addition, the use of piezoelectric elements for that matter grew in relevance, for its capability of creating those effects without penalizing other structural parameters, such as overall mass loading (Zhang et al. 2015;Nouh, Aldraihem, & Baz 2016;Zientek, Gardonio, & Dal Bo 2018).
When dealing with subwavelength resonant type metamaterials, the structural design of each cell can be intricate which, together with the necessity of manufacturing enough cells for assembling the full lattice structure, often leads to additive manufacturing (Bilal, Foehr, & Daraio 2017;Qureshi, Li, & Tan 2016).These techniques are, on their own, also gaining interest in many engineering areas, from material sciences and manufacturing, to those related to the performance assessment of such 3D printed parts.Among others, the Fused Deposition Modeling (FDM) is interesting for its capability of creating somehow complex geometries at a relatively low cost (Schumacher et al. 2015, Bosqué 2015).
Therefore, if metamaterials are to be incorporated into new products development cycles, there has to be reliable methods for design and simulation that can predict full and subsystems dynamic behavior and performance.In this way, as far as the design of FDM manufactured metamaterials is concerned, the accurate modeling of such structures is one of the key steps on the development of a metamaterial (sub)structures.The aim of using FDM parts together with smart materials, towards smart metamaterials, poses two main challenges: (i) the parts created by FDM processes may exhibit anisotropic behavior (Lovo, Fortulan & da Silva, 2018) and (ii) the electromechanical coupling, to which analytical methods are available through rather restrictive in term of geometry complexity.This paper addresses these issues, presenting the theoretical and practical aspects of the metamaterial modeling in section 2, followed by the experimental validation in section 3 and, finally, drawing some conclusions in section 4.

METASTRUCTURE MODELING
This section deals with the procedure adopted for the full system modeling.The first steps involve the modeling of a single unit cell, starting from mechanical and electrical subsystems and their integration into an electromechanical system.The state-space representation of the unit cells allows the assembly of a periodic arrangement of unit cells into the host structure, referred to as metastructure.
The unit cell concept is often used for the metamaterial analysis (Raghavan & Phani 2013;Claeys et al. 2015;Li, Zhang, & Liu 2017), as the full lattice is often composed by a periodic arrangement of these cells.Therefore, the modeling of these cells is a key point for the analysis of the metastructure.For an initial assessment of the cell's dynamic behavior, it is very common to use lumped parameter models, either purely mechanical or electromechanical as shown in Figure 1 (Hu et al., 2016;Sugino et al., 2016).One can write for a flexural electromechanical beam structure where m , c and k are the structural equivalent mass, damping and stiffness, respectively, x is the generalized displacement, v is voltage across the piezoelectric element,  is the electromechanical coupling and f is the external force acting on the unit cell.Equation ( 1) is often called the sensor equation and is related to the actuator equation (2), which describes the electric behavior of the system, with respect to the load circuit impedance Z and piezoelectric capacitance p C , Modular modeling approach for FDM printed structures and piezo disks for metamaterial design Gabriel K. Rodrigues et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e204 3/14 The electromechanical coupling becomes evident with the term dependent on x also present in Eq. ( 2).The piezoelectric system needs to be connected to a load circuitry, such as a LR shunt, which will provide the impedance Z and couple both domains.

Mechanical subsystem modeling
A lumped model can be used to describe the mechanical behavior of the system and would, arguably, be the most compact representation for such a system.Considering a homogeneous Euler-Bernoulli beam, one can write the potential energy U for the continuous system as: where the term EI stands for the flexural stiffness of the Euler-Bernoulli beam, L is the total length and the lateral displacement function of space x and time t is   , w x t .Also, the kinetic energy can be obtained with, in which  is the density of the beam and S is the cross section area.
As stated initially, a lumped model will be used, hence the finite form of Eqs. ( 3)-( 4) with n elements yields: Therefore, the stiffness and mass matrices can be written with Eqs. ( 5)-( 6) to achieve the elementary parameters for a discrete beam, as proposed by Eddanguir and Benamar (2013).As an alternative, the Finite Element Method (FEM) is often used for more complex problems, with arbitrary geometries and complex material behavior.One can use the interpolation provided by the FEM and write, with / l L n  being the element length, h the cross section height and b its width (c.f. Figure 2a), the elementary mass matrix e m and elementary stiffness matrix e k (Craig, 1981)    Both, analytical and finite element models, will be compared in the results section of this work.Either way, in this work, the mechanical state-space matrix is formulated in the modal domain i.e: where the I is the identity matrix 2 m W , is the diagonal matrix of squared natural frequencies and m D is the diagonal modal damping matrix for the mechanical system.
As the proposed material for the application is a FDM polymer, it is important to discuss the hypothesis used on the modeling process.As mentioned previously, a FDM produces an anisotropic structure (Lovo, Fortulan & da Silva, 2018).Even though, it is reasonable to assume an orthotropic material with planar isotropy, as the resistance in the printing plane does not vary significantly as compared with the resistance between layers.The porosity parameter can be taken into account by considering a penalty on the stiffness and/or density of the material.As this factor may vary between printers, printing materials and printing parameters, one should always consider calibrating model parameters against simple models or update these parameters based on preliminary prototypes, which is the case in this paper.

Electrical subsystem modeling
The electrical behavior of piezoelectric materials can be described as an electrical source connected to a capacitor, c.f. Figure 1(b).The external load impedance ( Z ) is often an inductance, a resistance or both.
As shown in Eq. ( 2), the first term indicates the current, i , through the load.Therefore, one can write for a pure resistance R , Similarly, for a pure inductance, J For a parallel association of a resistor and an inductor, a combination of Eqs. ( 10) and (11) yields: Substituting Eq. ( 12) in ( 2) and rearranging to a differential equation notation results in Modular modeling approach for FDM printed structures and piezo disks for metamaterial design Gabriel K. Rodrigues et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e204 5/14 or, The first term on the right hand side of Eq. ( 14) is the resonance frequency of the electrical system, the second one is related with damping and the third is the electromechanical coupling.Therefore, the electrical state-space matrix can be written as where W is the squared resonance frequency of the electrical model and the e D is the damping.One can easily read that if there is no resistor or inductor on the electrical circuit, the term vanishes from Eq. ( 13).

Electromechanical coupling
As described by Erturk & Inman (2008) as an analytic solution for a unimorph beam, the bending moment due to the piezoelectric patch is written as The total moment of a generic piezo with arbitrary length and position attached to a structure is given by Eq. ( 19), which includes the multiplication by Heaviside functions H : Therefore, one can fit Eq. ( 19) into a vector as: Modular modeling approach for FDM printed structures and piezo disks for metamaterial design Gabriel K. Rodrigues et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e204 6/14 where Q is the electromechanical coupling vector.Converting Eqs. ( 1)-( 2) into matrix form and arranging for state-space representation yields: Then it is possible to build the full state matrix of the unit cell as:

Full metastructure modeling
For the sake of simplicity, a flexural beam will be adopted as host structure, with structural state space matrix in which, W is the diagonal matrix of squared natural frequencies and host D is the diagonal modal damping matrix of the host structure.The strength of this technique resides exactly in having all elements of the metastructure chain represented in the same fashion.Each state-space modal model can be derived from different methodologies, depending on the governing physical phenomena and the more convenient modeling tool at hand, such as analytical models, lumped parameter, finite elements, experimental modal analysis, etc.
As the cells and the host subsystems will share some degree of freedom, thanks to the mechanical attachments amongst each other, the off diagonal matrixes ch,i A enables the links between these modules, resulting in the full metastructure state-space dynamic matrix in (25).

Experimental setup
For the experimental validation, the selected benchmark structure consists of a FDM beam equipped with a piezo disk (Figure 3), that way, both analytical and FE models can be used for validation purposes.Base excitation, with an electrodynamic shaker (model K2007E01 from Modal Shop), was used as input, while the outputs were the tip velocity and the voltage across the piezo electrodes.The force input was measured with a 208C02 force transductor from PCB Piezotronics; the tip velocity was measured by a Laser Doppler Vibrometer CLV1000 by Polytec.All the sensors, including the voltage measurements, were acquired with the LMS SCADAS mobile system running LMS TestLab for processing the data, then exported to MATLAB for post-processing and data plotting.On the electrical connections, a resistive and an inductive circuit were used (Figure 4).The inductive circuit was tuned to match the first mechanical mode with the use of a synthetic inductor.A synthetic inductance was preferred on this test, as the inductance for low frequency resonance on an electrical system is reasonably high (above 100F), as shown in Table 1.

Experimental results
Figure 5 shows structural frequency response functions (FRFs) from base input force to the beam's tip velocity.Four circuit configurations are tested, namely, only the piezo in open and short circuit, resistive and inductive impedances.As shown, there is not a significant reduction on tip velocity due to circuitry manipulation.This suggests a rather weak coupling between the FDM beam and the piezo elements adopted.Even though, it is noticeable that the shunt circuit is performing as expected, and a minor amplitude reduction is present.Figure 6 shows the FRFs from force input to voltage output (c.f.0 v in Figure 4).The difference in electrical response between a tuned and a detuned shunt circuit is clear.While the detuned system shows a resonance around 9Hz and a rather flat response over the frequency range, the tuned circuit shows a drastic difference between resonant and offresonance behavior, with a five fold increase in the peak value close to 100Hz.This indicates that if the selected piezo patch presented a stronger coupling with the host structure, a reasonable performance in mechanical peak suppression would have been present.Also, a pure resistive load was compared against the tuned inductive load, to which results are depicted in Figure 7.Even thou the resistive load used is not the optimum value for this situation, its peak value is two orders of magnitude below the tuned one.Further discussion on this topic will follow in the next Section.

Model validation
The electromechanical models, as described in Section 2, have been implemented in MATLAB.In order to allow a direct comparison, the experimental data also comprises the chosen outputs as the tip velocity and electric voltage.
Figure 8 shows the tip velocity from both lumped and FEM models, in order to crosscheck both modeling strategies.This comparison is based on driving-point FRFs, when excitation and response are taken in the same degree of freedom, in this case the beam tip.The comparison shows that both models deliver very similar values of amplitude throughout the frequency range.
Similarly, Figure 9 presents results for the electrical output, when a resistive load is used.As it can be seen both models agree.The higher electrical response is expected, since the structural amplitude on the FEM model is also slightly higher.From this point on, therefore, the FEM model with base excitation is used as a reference for the experimental validation.
While Figure 10 shows a comparison for the beam tip velocity, Figure 11 shows the electrical response across the piezoelectric element.It is possible to conclude that the modeling procedure is effective by comparing the FEM model to the experimental data.Also in Figure 11, the maximum theoretical electrical output for a resistive load is shown.Its peak value, around 6V/N, sits closely below the detuned inductor, around 10V/N, but stills fairly below the tuned inductor, around 50V/N.
The dynamics observed on the experimental data, below the predicted 97Hz resonance frequency, is due to the fixture, which is not an ideal clamping mechanism, as assumed in the FEM model.The offset between the experimental and numerical resonance frequency is attributed both, to the clamping as well as to the effect of porosity, neglected on the modeling process.Even thou, the results come quite close together.

Metastructure simulation
As a proof of concept a full metastructure was simulated with the proposed method, using the models previously validated for the smart unit cell.The full system consists of a beam host structure, clamped in both ends and equipped with periodically placed beam like resonators, see Figure 12.The overall added mass to the host structure is around 21%.Both input (force) and output (velocity) are situated on the middle point of the host beam.As shown in Figure 13, the mechanical behavior of the resonators are tuned for the first mode of the host structure, thus, they create a band gap around that frequency.In addition, the electric behavior of the resonant circuit was tuned to tackle the border peaks that emerge from band the gap creation, rather than the bandgap itself.Thus, they are responsible for the attenuation of the vicinity peaks and, in practice would widen the range in which the overall vibration of the host structure is reduced.
The electrical resonators were interleaved leading to a "ABAB" pattern, where an "A" cell has electrical tuning for the lower frequency peak (before the band gap), while the "B" cell is electrically tuned for the higher peak (after the band gap).The mechanical behavior of such a single cell is displayed on Figure 13, with tip mobility (due to a tip force input) and a short circuited piezo disk.In this way, only the mechanical effect of the piezo is present.
A 35Hz bandgap is generated around the structure resonance, however, rather tall sideband peaks results.As it can be seen in Figure 13, this novel strategy of targeting sidepeaks with smart unit cells is promising, as a clear attenuation is present, even considering piezoelectric elements with a low coupling.The combined effects of the number of cells present on the periodic structure allows for it, resulting in some 30dB reduction on structural vibration on the borders of the bandgap (dashed line in Figure 13).

CONCLUSIONS
A modular modeling approach for electromechanical systems has been presented in this paper.As such, it offers an accurate and versatile modeling strategy for smart unit cells built via FDM process.The model is capable of using both lumped and FEM strategies with ease, without losing its accuracy.Arguably, it even allows for hybrid simulation, when subsystem models are based on experimental data (such as modal analysis).
The use of the proposed modular approach, based on a state-space formulation, facilitates the coupling of multiphysics elements on the unit cell level, namely, the structural resonator, the piezoelectric element and load circuit.The electromechanical coupling between the two domains is retained in the state-space form, which is key for both direct and indirect piezoelectric effect exploitation.
The modal state-space formulation of the unit cells also allows some independency between the modeling of the host structure and the electromechanic unit-cell, as well as the assembly of the full metastructure.This method also enables the use of different types of cells with different levels of discretization, e.g. a host structure with very coarse lumped model and very refined lumped resonator.
Besides validating the proposed approach, this work also presented case study, based on a clamped-clamped metallic host beam installed with FDM printed smart unit cells.A novel approach for exploiting low efficient piezoelectric elements is presented, showing promising results.It is based on an "ABAB" pattern for the periodic elements that target the side peaks that rise from the bandgap creation rather than the bandgap itself.Numerical results indicate a potential for vibration attenuation at those frequencies, which widens the useful range for the metastructure.

Figure 2 :
Figure 2: Element parameters for host structure and piezoelectric element (a) structural element (b) unimorph parameters and input moment.

Y
is the Young modulus of the piezo patch, p b the width of the piezo, p h the piezo height, 1 h and 2 h are respectively the distance from the bottom and top of the piezo to the host structure neutral axis, 31 d is a piezoelectric constant expressed in coulomb per newton [C/N], and the spatial integral over y corresponds to the electromechanical coupling  , which can be evaluated as the distance between the neutral axis of the piezo to the neutral axis of the beam, and 31 e is another piezoelectric constant in coulomb per square meter [C/m 2 ] correlated with 31 d by the Young modulus.

Figure 3 :
Figure 3: Test setup for model validation including FDM beam with piezo disk base driven by a shaker

Figure 5 :
Figure 5: Mechanical response comparison for different circuit configurations

Figure 7 :
Figure 7: Electrical response -influence of load nature

Figure 12 :
Figure 12: Proposed metastructure with periodically beam like resonators

Figure 13 :
Figure 13: Mobility of structures -simple resonator to full metastructure simulation as:

Table 1 :
Electrical system parameters