Investigation and implementation of Active Willis Metamaterials (AWM) have been done exclusively, in all the available literature, by approaches that do not rely on any solid control theory basis. When coupled with piezoelectric control elements, the available approaches have not included, from the first principles, the exact form of the constitutive relationship of the piezoelectric materials. Furthermore, in all these approaches, stability analysis, robustness, ability to accommodate uncertainty or parameter changes, or consideration of disturbance rejection has not been addressed at all. More importantly, the available formulations have always mixed the flow and effort variables of the AWM, resulting in a form that is totally incompatible for the use in generating, investigating, or even designing any appropriate sensing or control applications of the material. In this paper, the piezoelectric-based AWM is modeled, from the first principles, to develop a constitutive coupling form that enables its use in actuation, sensing, and as an integrated controller that can be analyzed, designed, and optimized using the classical, optimal, and robust control system theories. Lagrange dynamics formulation is used to generate the equations governing the Willis coupling, the piezoelectric coupling, as well as the active robust controller. With this developed controlled-based structure of the AWM, the inherent and powerful capabilities of the AWM that lie in its ability to robustly control the material properties themselves such as the compliance (or stiffness) and specific volume (or density) are demonstrated in great detail via several numerical examples. Controlling these properties enables the AWM to be used in numerous important and imaginative applications such as cloaking, beam shifting, beam focusing, as well as many other applications that are limited only by our imagination.

Willis materials, since their inception by Willis,1 have been the focus of extensive research, particularly during the past few years because of their unique and inherent capabilities in both their passive and active forms. In their passive form, these materials exhibit elasto-dynamic behavior that couples the stress and momentum to the strain and particle velocity. This unique coupling characteristic arises from both local effects and material asymmetries.2–4 For example, Muhlestein2 and Muhlestein et al.3 demonstrated experimentally the Willis material behavior by using an air duct with asymmetric restriction that is membrane-loaded from one direction. Also, Melnikov et al.4 developed an innovative acoustic meta-atom experiment to verify the existence of the maximum Willis coupling phenomenon. Furthermore, Groby et al.5 investigated the Willis material behavior and extracted its parameters using air ducts coupled with detuned Helmholtz resonators, which are located at different positions along the ducts. Distinct also among the experimental validation of the Willis material behavior is the work of Liu et al.,6 whereby beams are provided with an array of asymmetrically designed discs to control the wave propagation along the forward and backward directions.

The theory that governs, explains the fundamentals, and sets the physical bounds on the attractive features of the above passive Willis materials’ demonstrations has been reported extensively in many accounts such as, for example, in the work of Milton and Willis,7 Pernas-Salomón and Shmuel,8 Sieck et al.,9 and Muhlestein et al.10,11

With such unique capabilities as those of the Willis materials, it has been shown that this class of materials can be employed in a wide variety of applications such as the design and of Parity-Time (PT) symmetric systems,12 shaping and manipulation of the acoustic wave front,13 as well as synthesizing non-reciprocal systems.14,15

It is important here to note that passive Willis materials, as mentioned above, are synthesized metamaterials where the coupling effects do not appear in the constitutive equations of the individual ingredients. But rather, the asymmetric integration of these ingredients generates the unique coupling effects that provide the Willis composite with its attractive behavior, which is completely different from that of the individual gradients. Such distinction can further be enhanced by providing the Willis composite with active control capabilities. Extensive efforts have been exerted along this direction,16–25 whereby the performance of the Willis material is enhanced and its behavior is tuned in order to adapt to the needs of the considered applications. In all these studies, the Willis material is integrated with piezoelectric materials that have natural built-in coupling characteristics that couple the mechanical stress and strain with the electrical field and electrical displacement. With the resulting Willis-Piezoelectric composite, it becomes possible to further couple the electrical field with the velocity and momentum fields, resulting in a unique electro-momentum coupling that would enable the sensing, actuation, and control of the composite in a way that creates a new class of tunable metamaterials’ properties.

In all the cited active Willis metamaterials (AWM), the constitutive equations are always written as:
{ σ D p } = [ C B S B A W S W ρ ] { ε E v } ,
where σ , D , a n d p a r e the stress, electric displacement, and momentum, respectively. Also, ε , E , a n d v a r e the strain, electric field, and velocity, respectively. Note that C and ρ denote the stiffness modulus and density, whereas B and W are the coupling coefficients between the piezo and Willis materials, which have the Willis coupling S. Also, A denotes the strain-free permittivity of the piezoelectric material.

Note that the left-hand-side (LHS) of the constitutive equation is a mix between effort variables (σ, p) and flow variable (D), whereas the right-hand-side (RHS) has also a mix between an effort variable (E) and flow variables (ɛ, v). This forms the main obstacle in using this material for active control applications. Furthermore, the coupling terms B and W are extracted by not considering the basic and fundamental form of the piezoelectric constitutive equations.

With such a developed control-based structure of the AWM, it is envisioned that this class of metamaterial will be able to have tunable material properties such as compliance (or stiffness) and specific volume (or density) as will be demonstrated in the present paper. Furthermore, it is also expected that the AWM will have controllable bandgap and dispersion characteristics that will enable its use in numerous important and imaginative applications.

Therefore, this paper aims at generating the constitutive equations of the AWM in a form that enables its use in actuation, sensing, and as an integrated controller that can be analyzed, designed, and optimized using the classical, optimal, and robust control system theories. Lagrange dynamics formulation is used to generate the equations governing the Willis coupling, the piezoelectric coupling, as well as the active robust controller.

Accordingly, the paper is organized in eight sections. In Sec. I, a brief introduction is presented. The concept of the active piezoelectric Willis metamaterial is introduced in Sec. II. The modeling of the active piezoelectric-Willis material is developed in Sec. III. The operation of the active Willis material (AWM) for sensing is presented in Sec. IV. The operation of the active Willis material (AWM) for control of flow variables is presented in Sec. V. Section VI presents the active Willis material (AWM) with controllable properties, and Sec. VII demonstrates the effectiveness of the AWM in generating controllable specific volume and compliance while accommodating parameter uncertainty and rejecting the effect of external disturbances. Section VIII summarizes the conclusions of the present study and the potential for its future extensions.

In this paper, the emphasis is placed on the active piezoelectric-Willis metamaterial configuration shown in Fig. 1 which has been analyzed extensively by Muhafra et al.18 In that configuration, two different piezoelectric linear springs are used to connect three different point masses to form an AWM unit cell. In this manner, this configuration presents a general arrangement whereby both the stiffness and the mass are asymmetrically distributed over the unit cell. Such a configuration has been inspired from the work of Muhlestein et al.3 who considered a much simpler configuration consisting of a linear spring linking two different point masses to investigate the reciprocity, passivity, and causality properties of the passive Willis materials.

FIG. 1.

Active periodic Willis material (AWM).

FIG. 1.

Active periodic Willis material (AWM).

Close modal

The considered AWM has been selected as it has been shown by Muhafra et al.18 to exhibit, on one hand, a coupling between the stress–velocity and the momentum–strain constitutive equations, which is known as the “Willis Coupling.” On the other hand, it has also been shown that there is a coupling between the electric displacement—velocity—and the electric field—momentum, which is known as the “Piezoelectric Coupling,” which is analogous to the Willis passive material coupling.

It is important to note that the nature of this generated Willis coupling is developed primarily from the material asymmetries.

With a material such as the AWM with built-in Willis and piezoelectric couplings, it is possible to utilize the resulting electro-momentum coupling as an effective and powerful degree of freedom for imparting sensing and actuation capabilities. More importantly, this additional degree of freedom can further enable the control and tuning of the physical properties of the AWM as needed for different applications.

The equations of motion of the system are generated by using the Lagrange's dynamics approach and the constitutive equations of the piezo-actuators, which have open-circuit stiffness k 1 Dand k 2 D.

The kinetic energy
K E = 1 2 m 1 u ¨ L 2 + 1 2 m 2 u ¨ m 2 + 1 2 m 3 u ¨ R 2 ,
(1)
where m1,..,3 are the masses of the three blocks shown in Fig. 1 and uL,m,R denote the translation of the left, middle, and right masses, respectively.
The potential energy (Details are given in  Appendix)
P E = 1 2 k 1 D ( u L u m ) 2 + 1 2 1 C p 1 s Q 1 2 h 1 Q 1 ( u L u m ) + 1 2 k 2 D ( u m u R ) 2 + 1 2 1 C p 2 s Q 2 2 h 2 Q 2 ( u m u R ) ,
(2)
where C p 1 sand C p 2 sare the strain-free capacitances of the two piezoelectric springs. Also, h1 and h2 are the corresponding piezoelectric coupling coefficients. Furthermore, Q1 and Q2 denote the piezoelectric charges associated with the two springs.
Hence, the Lagrangian L is given by
L = K E P E .
(3)

Also, the generalized forces are given by FL and FR acting on generalized DOF uL and uR and also V1 and V2 acting on generalized DOF Q1 and Q2.

Accordingly, the Equations of Motion are extracted using Lagrange's dynamics as follows:

Structural:
d d t ( L u ˙ L ) L u L = F L , d d t ( L u ˙ m ) L u m = 0 , d d t ( L u ˙ R ) L u R = F R .
Piezoelectrical:
d d t ( L Q ˙ 1 ) L Q 1 = V 1 and d d t ( L Q ˙ 2 ) L Q 2 = V 2 .

This leads to,

Structural:
m 1 u ¨ L + k 1 D ( u L u m ) h 1 Q 1 = F L ,
(4)
m 2 u ¨ m k 1 D ( u L u m ) k 2 D ( u R u m ) + h 1 Q 1 h 2 Q 2 = 0 ,
(5)
m 3 u ¨ R + k 2 D ( u R u m ) + h 2 Q 2 = F R .
(6)
Piezoelectrical:
1 C p 1 s Q 1 h 1 ( u L u m ) = V 1
(7)
and
1 C p 2 s Q 2 + h 2 ( u R u m ) = V 2 .
(8)
Equations (7) and (8) give
h 1 Q 1 = h 1 C p 1 s V 1 + C p 1 s h 1 2 ( u L u m )
(9)
and
h 2 Q 2 = h 2 C p 2 s V 2 C p 2 s h 2 2 ( u R u m ) .
(10)
Equations (5), (9), and (10) give
( m 2 s 2 + k 1 D + k 2 D C p 1 s h 1 2 C p 2 s h 2 2 ) u m = ( k 1 D C p 1 s h 1 2 ) u L + ( k 2 D C p 2 s h 2 2 ) u R C p 1 s h 1 V 1 + C p 2 s h 2 V 2
or
u m = [ ( k 1 D C p 1 s h 1 2 ) u L + ( k 2 D C p 2 s h 2 2 ) u R C p 1 s h 1 V 1 + C p 2 s h 2 V 2 ] ( m 2 s 2 + k 1 D + k 2 D C p 1 s h 1 2 C p 2 s h 2 2 ) .
In the dimensionless form,
u m = [ ( 1 κ 1 2 ) u L + k r ( 1 κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] [ μ 2 s ¯ 2 + ( 1 κ 1 2 ) + k r ( 1 κ 2 2 ) ] ,
(11)
where κ 1 2and κ 2 2are the electromechanical coupling factors of the first and second piezoelectric actuators. Also, s = Laplace complex number, k r = k 2 D / k 1 D μ2 = m2/m1, ω n 2 = k 1 D / m 1, and s ¯ = s / ω n. Furthermore, d1 and d2 denote the piezoelectric strain coefficients of the first and second piezoelectric actuators.
Substituting Eqs. (9)–(11) into Eq. (4) gives
( m 1 s 2 + k 1 D C 1 s h 1 2 ) u L ( k 1 D C 1 s h 1 2 ) [ ( 1 κ 1 2 ) u L + k r ( 1 κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] [ μ 2 s ¯ 2 + ( 1 κ 1 2 ) + k r ( 1 κ 2 2 ) ] h 1 C p 1 s V 1 = F L .
In the dimensionless form,
( s ¯ 2 + ( 1 κ 1 2 ) ) u L ( 1 κ 1 2 ) [ ( 1 κ 1 2 ) u L + k r ( 1 κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] [ μ 2 s ¯ 2 + ( 1 κ 1 2 ) + k r ( 1 κ 2 2 ) ] ( 1 κ 1 2 ) d 1 V 1 = 1 k 1 D F L
or
( s ¯ 2 + ( 1 κ 1 2 ) ( 1 κ 1 2 ) 2 D ) u L k r ( 1 κ 1 2 ) ( 1 κ 2 2 ) D u R + ( ( 1 κ 1 2 ) 2 D + ( 1 κ 1 2 ) ) d 1 V 1 k r ( 1 κ 1 2 ) ( 1 κ 2 2 ) D d 2 V 2 = 1 k 1 D F L ,
(12)
where D = [ μ 2 s ¯ 2 + ( 1 κ 1 2 ) + k r ( 1 κ 2 2 ) ].
Substituting Eqs. (9) through (11) into Eq. (6) gives
( m 3 s 2 + k 2 D + C p 2 s h 2 2 ) u R ( k 2 D C p 2 s h 2 2 ) [ ( 1 κ 1 2 ) u L + k r ( 1 + κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] D + h 2 C p 2 s V 2 = F R .
In the dimensionless form,
( μ 3 s ¯ 2 + k r ( 1 + κ 2 2 ) k r 2 ( 1 κ 2 4 ) D ) u R k r ( 1 κ 1 2 ) ( 1 κ 2 2 ) D u L + k r ( 1 κ 1 2 ) ( 1 κ 2 2 ) D d 1 V 1 ( k r 2 ( 1 κ 2 2 ) 2 D k r ( 1 κ 2 2 ) ) d 2 V 2 = 1 k 1 D F R .
(13)

In a matrix form,

(14)

In a compact form, Eq. (14) reduces to

(15)
The transformation from force
{ F L F R } T
 - deflection { u L u R } Tdomain to the local stress-momentum { σ p } T - strain-velocity { ε υ } T domain is carried out such that
(16)
where L is the distance between any two masses at equilibrium.
Substituting Eq. (16) into Eq. (15) gives
[ K D ] [ M ( 2 ) ] 1 { ε υ } = 1 k 1 D [ M ( 1 ) ] 1 { σ p } + [ H ] { d 1 V 1 d 2 V 2 }
(17)
or
{ ε υ } = 1 k 1 D ( [ M ( 2 ) ] [ K D ] 1 [ M ( 1 ) ] 1 ) { σ p } + ( [ M ( 2 ) ] [ K D ] 1 [ H ] ) { d 1 V 1 d 2 V 2 } = s V { σ p } + d { d 1 V 1 d 2 V 2 } ,
(18)
where s V = 1 k 1 D ( [ M ( 2 ) ] [ K D ] 1 [ M ( 1 ) ] 1 ) = short-circuit compliance and density matrix and d = ( [ M ( 2 ) ] [ K D ] 1 [ H ] ) = piezo-strain coefficients matrix.
Equations (9) and (10) are now modified such that
Q 1 = C p 1 s h 1 u L C p 1 s h 1 [ ( 1 κ 1 2 ) u L + k r ( 1 κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] D + C p 1 s V 1 = C p 1 s h 1 ( 1 ( 1 κ 1 2 ) D ) u L C p 1 s h 1 ( k r ( 1 κ 2 2 ) D ) u R + ( C p 1 s h 1 ( 1 κ 1 2 ) D C p 1 s d 1 ) d 1 V 1 C p 1 s h 1 k r ( 1 κ 2 2 ) D d 2 V 2
(19)
and
Q 2 = C p 2 s h 2 u R + C p 2 s h 2 [ ( 1 κ 1 2 ) u L + k r ( 1 κ 2 2 ) u R ( 1 κ 1 2 ) d 1 V 1 + k r ( 1 κ 2 2 ) d 2 V 2 ] D + C p 2 s V 2 = C p 2 s h 2 ( 1 κ 1 2 ) D u L C p 2 s h 2 ( 1 k r ( 1 κ 2 2 ) D ) u R C p 2 s h 2 ( 1 κ 1 2 ) D d 1 V 1 + ( C p 2 s h 2 k r ( 1 κ 2 2 ) D + C p 2 s d 2 ) d 2 V 2 .
(20)

Combining Eqs. (19) and (20) in a matrix form gives

(21)
Equation (21) is rewritten as
(22)
Equation (22) can be rearranged as follows:
(23)
where = charge sensitivity matrixand = piezo-capacitance matrix.
Equations (18) and (23) are combined and arranged in the following matrix form:
(24)

Then, the final form of the Willis Material model is given by

(25)

Equation (25) indicates that the LHS

{ ε υ Q 1 Q 2 } T
defines the output vector of the structural flow variables { ε υ } T followed by the electrical flow variables { Q 1 Q 2 } T. On the other hand, the RHS { σ p d 1 V 1 d 2 V 2 } T defines the input vector of the structural effort variables { σ p } T followed by the electrical effort variables { d 1 V 1 d 2 V 2 } T. This particular formulation of the Active Willis Material (AWM) is particularly suitable for systematically structuring active control systems that couple further the input effort variables to the output flow variables via appropriate feedback strategies.

Figure 2 displays a visualization of the inherent features that govern the interactions between the effort and flow variables of the Active Willis Material (AWM) as described in Eq. (25).

FIG. 2.

The inherent coupling structure of the active Willis material (AWM).

FIG. 2.

The inherent coupling structure of the active Willis material (AWM).

Close modal

Specifically, Fig. 2 displays, in its LHS, the nature of interaction between the structural and electrical effort variables { σ p d 1 V 1 d 2 V 2 } T on one hand and the structural flow variables { ε υ } T on the other hand. Such an interaction is described in the top half of Eq. (25).

Furthermore, it displays, in its RHS, the corresponding interaction with the electrical flow variables { Q 1 Q 2 } T as described in the bottom half of Eq. (25).

Figure 3 casts this interaction in the form of the open-loop configuration of the Active Willis Material (AWM).

FIG. 3.

Open-loop configuration of the AWM.

FIG. 3.

Open-loop configuration of the AWM.

Close modal

This mode of operation requires that the piezo-springs operate under short-circuit conditions.26–28 

  • Sensing of {σ p}

The bottom half of Eq. (25) reduces to
{ σ p } = [ g s V ] 1 { Q 1 Q 2 } .
(26)

Equation (26) is represented graphically by the block diagram of Fig. 4(a).

  • Sensing of {ɛ υ}

FIG. 4.

The active Willis material (AWM) as a sensor. (a) Sensing {σ p}. (b) Sensing {ɛ υ}.

FIG. 4.

The active Willis material (AWM) as a sensor. (a) Sensing {σ p}. (b) Sensing {ɛ υ}.

Close modal
Combining the top half of Eq. (25) and Eq. (26) gives
{ ε v } = [ s V d ] { σ p } = [ s V d ] [ g s V ] 1 { Q 1 Q 2 } .
(27)

Equation (28) is represented graphically by the block diagram of Fig. 4(b).

In this mode of operation, the AWM is used to control the structural flow variables {ɛ v} by generating the control voltages V1 and V2 according to the following control law:
{ d 1 V 1 d 2 V 2 } = [ G 1 ε G 1 v G 2 ε G 2 v ] ( { ε R v R } { ε v } ) = [ G ] ( { ε R v R } { ε v } ) ,
(28)
where [G] is the control gain matrix and { ε R υ R } T is a reference and desired set of the structural flow variables {ɛ υ}.
Using the top half of Eq. (25) along with Eq. (28) yields
{ ε v } = s V { σ p } + d [ G ] ( { ε R v R } { ε v } )
or
{ ε v } = ( I + d [ G ] ) 1 ( s V { σ p } + d [ G ] { ε R v R } ) .
(29)

Equation (29) is illustrated graphically by the block diagram of Fig. 5.

FIG. 5.

Closed-loop configuration of the AWM.

FIG. 5.

Closed-loop configuration of the AWM.

Close modal

Note that the control gain matrix [G] can be determined using the robust control theory to reject the effects of the external disturbances, which are generated by the effort variables {σ p}.26–28 

Such an application is rather simple and does not demonstrate the inherent and powerful capabilities of the AWM, which lie in its ability to control the material properties themselves such as the stiffness and density. Controlling these properties enables the AWM to be used in numerous important and imaginative applications such as cloaking, beam shifting, beam focusing, as well as many other applications that are limited by our imagination. The following section will present the theoretical foundation for controlling these properties and demonstrate, by numerical examples, the potential behavior of AWM provided with such controllable and programmable characteristics.

In this section, the density and stiffness of the AWM are controlled.

The top half of Eq. (25) is expanded as follows:
{ ε v } = s V { σ p } + d [ d 1 V 1 d 2 V 2 ] = [ s 11 V s 12 V s 21 V s 22 V ] { σ p } + [ d 11 d 12 d 21 d 22 ] [ d 1 V 1 d 2 V 2 ] .
(30)
The second row of Eq. (30) is expanded to give
v = [ s 21 V s 22 V ] { σ p } + [ d 21 d 22 ] [ d 1 V 1 d 2 V 2 ] = s 21 V σ + s 22 V p + [ d 21 d 22 ] [ d 1 V 1 d 2 V 2 ] .
(31)
Dividing by p gives the closed-loop specific volume υ (i.e., the reciprocal of the density) as follows:
υ = v p = 1 ρ C l o s e d L o o p S p e c i f i c V o l u m e = s 22 V O p e n L o o p S p e c i f i c V o l u m e + 1 p [ d 21 d 22 ] [ d 1 V 1 d 2 V 2 ] C o n t r o l l e r o f S p e c i f i c V o l u m e + s 21 V σ p E f f e c t o f E x t e r n a l D i s t u r b a n c e .
(32)

Note that the present control-based formulation of the Active Willis Material results in casting its properties in terms of the specific volume and compliance instead of the density and the modulus.

Equation (32) indicates that in order to obtain a desired specific volume υ R, it is a must to design a controller (d1V1 and d2V2) that can reject the effect of the external disturbance s 21 V ( σ / p ) and accommodate the uncertainty generated by including p in the controller term 1 p [ d 21 d 22 ] [ d 1 V 1 d 2 V 2 ] T.

If the control voltages, [ d 1 V 1 d 2 V 2 ] T are selected such that
[ d 1 V 1 d 2 V 2 ] = [ g 1 g 2 ] ( υ R υ ) = g υ [ 1 g r υ ] ( υ R υ ) ,
(33)
where g1 and g2 are the control gains for d 1 V 1 a n d d 2 V 2, respectively. Also, gυ = g1 and g r υ = g 2 g υ. Then, Eq. (32) can be rewritten as
( υ R υ ) = g υ p ¯ ( d 21 + g r υ d 22 ) ( υ R υ ) + ( υ R s 22 V ) s 21 V σ p .
(34)

In Eq. (34), p ¯ denotes a nominal value of the momentum to be estimated at the nominal operating frequency of the AWM.

Manipulation of Eq. (34) yields
( υ R υ ) T o t a l E r r o r = [ 1 1 + K υ ] υ R T r a c k i n g E r r o r + ( [ 1 1 + K υ ] s 22 v + [ 1 1 + K υ s 21 v ] σ p ) E r r o r d u e t o E x t e r n a l D i s t u r b a n c e ,
(35)
where K υ = g υ p ¯ ( d 21 + g r υ d 22 ).
Equation (35) can be written as
( υ R υ ) = [ G 1 υ υ R + ( G 1 υ s 22 v + G 2 υ σ p ) ] ,
(36)
where G 1 υ = [ 1 1 + K υ ] and G 2 υ = [ 1 1 + K υ s 21 v ].

In Eqs. (35) and (36), the components of the total error are identified as the tracking error and error due to external disturbances. The controller gain Kυ, with components gυ and g r υ, is selected ideally to minimize the tracking error and, at the same time, reject the effect of the external disturbances. This minimization requires that Kυ should be as large as implied by Eq. (35). However, this leads to an excessively large control effort K υ ( υ R υ ). Therefore, an effective design of the control gain should ensure a balance that minimizes simultaneously the sum of the errors and the control effort.

But as the AWM in its open-loop formulation is undamped, the gain Kυ is selected to be in the form of a Proportional-Integral-Derivative controller to add the necessary damping such that
K υ = k p υ ( 1 + k ¯ υ i / s + k ¯ υ d s ) ( d 21 + g r υ d 22 ) 1 p ¯ ,
(37)
where k p υ = P r o p o r t i o n a l g a i n , k υ i = k p υ k ¯ υ i = I n t e g r a l g a i n , a n d k υ d = k p υ k ¯ υ d = D e r i v a t i v e g a i n .
Accordingly, the control effort can be expressed as follows:
E υ = K υ ( υ R υ ) .
(38)
Substituting Eq. (36) into Eq. (38) gives
E υ = K υ [ G 1 υ υ R + ( G 1 υ s 22 v + G 2 υ σ p ) ] .
(39)
Hence, the problem of optimal selection of the gain K υ is formulated as follows:
Find k p υ , k ¯ υ i , k ¯ υ d , a n d g r υ To minimize H υ 2 2 + E υ 2 2 = 1 2 π [ 0 ( | G 1 υ ( i ω ) | 2 + | K υ ( i ω ) | 2 ) ( | υ R | 2 + | s 22 v | 2 + + | s 21 v σ p | 2 ) d ω ] such that the closed loop system is stable ,
(40)
where H υ 2 and E υ 2 denote the H-2 norms of the specific volume and associated control effort, respectively.
The structure of the stiffness controller is extracted in a similar manner to that of the density controller described above. In this case, the first row of Eq. (30) is expanded to give
ε = [ s 11 V s 12 V ] { σ p } + [ d 11 d 12 ] [ d 1 V 1 d 2 V 2 ] = s 11 V σ + s 12 V p + [ d 11 d 12 ] [ d 1 V 1 d 2 V 2 ] .
(41)
Dividing by σ gives the closed-loop compliance s (i.e., the reciprocal of the modulus) as follows:
s = ε σ C l o s e d L o o p C o m p l i a n c e = s 11 V O p e n L o o p C o m p l i a n c e + 1 σ [ d 11 d 12 ] [ d 1 V 1 d 2 V 2 ] C o n t r o l l e r o f C o m p l i a n c e + s 12 V p σ E f f e c t o f E x t e r n a l D i s t u r b a n c e .
(42)

Equation (40) indicates that in order to obtain a desired compliance s R, it is a must to design a controller (d1V1 and d2V2) that can reject the effect of the external disturbance s 12 V ( p / σ ) and accommodate the uncertainty generated by including σ in the controller term 1 σ [ d 11 d 12 ] [ d 1 V 1 d 2 V 2 ] T.

If the control voltages [ d 1 V 1 d 2 V 2 ] T are selected such that
[ d 1 V 1 d 2 V 2 ] = [ g 3 g 4 ] ( s R s ) = g s [ 1 g r s ] ( s R s ) ,
(43)
where g3 and g4 are the control gains for d 1 V 1 a n d d 2 V 2, respectively, also, gs = g3 and g r s = g 4 / g 3 ., then, Eq. (43) can be rewritten as
( s R s ) = g s σ ¯ ( d 11 + g r s d 12 ) ( s R s ) + ( s R s 11 V ) s 12 V p σ
(44)
or
( s R s ) T o t a l E r r o r = 1 1 + K s ( s R s 11 v ) T r a c k i n g E r r o r + 1 1 + K s s 12 v p σ E r r o r d u e t o E x t e r n a l D i s t u r b a n c e ,
(45)
where K s = g s σ ¯ ( d 11 + g r s d 12 ).
Equation (45) can be rewritten as
( s R s ) = G 1 s ( s R s 11 v ) + G 2 s p σ ,
(46)
where G 1 s = [ 1 1 + K s ] and G 2 s = 1 1 + K s s 12 v.

In Eqs. (45) and (66), the components of the total error are identified as the tracking error and error due to external disturbances. The controller gain Ks, with components gs and g r s, is selected ideally to minimize the tracking error and, at the same time, reject the effect of the external disturbances. This minimization requires that Ks should be as large as implied by Eq. (45). However, this leads to an excessively large control effort K s ( s R s ). Therefore, an effective design of the control gain should ensure a balance that minimizes simultaneously the sum of the errors and the control effort.

But as the AWM in its open-loop formulation is undamped, the gain Ks is selected to be in the form of a Proportional-Integral-Derivative controller to add the necessary damping such that:
K s = k p s ( 1 + k ¯ s i / s + k ¯ s d s ) ( d 11 + g r s d 12 ) 1 σ ¯ ,
(47)
where k p s = P r o p o r t i o n a l g a i n , k s i = k p s k ¯ s i = I n t e g r a l g a i n , a n d k s d = k p s k ¯ s d = D e r i v a t i v e g a i n .
Accordingly, the control effort can be expressed as follows:
E s = K s ( s R s ) .
(48)
Substituting Eq. (46) into Eq. (48) gives
E s = K s [ G 1 s s R + ( G 1 s s 11 v + G 2 s p σ ) ] .
(49)
Hence, the problem of optimal selection of the gain K s is formulated as follows:
Find k p υ , k ¯ υ i , k ¯ υ d , a n d g r υ To minimize H s 2 2 + E s 2 2 = 1 2 π [ 0 ( | G 1 s ( i ω ) | 2 + | K s ( i ω ) | 2 ) ( | s R | 2 + | s 11 v | 2 + + | s 12 v p σ | 2 ) d ω ] such that the closed loop system is stable ,
(50)
where H s 2 and E s 2 denote the H-2 norms of the compliance and associated control effort, respectively.

Consider a Willis material (WM) with the physical properties listed in Table I.

TABLE I.

Physical properties of the Willis material.

Property k 1 D (n/m)kr C p 1 s (μf) C p 2 s (μf)μ2μ3ωn1 (rad/s) κ 1 2 κ 2 2d1 (m/volt)d2 (m/volt)L (m)
Value 1000 0.25 1.5 0.7 0.9 100 0.4 0.4 120 × 10−12 120 × 10−12 0.1 
Property k 1 D (n/m)kr C p 1 s (μf) C p 2 s (μf)μ2μ3ωn1 (rad/s) κ 1 2 κ 2 2d1 (m/volt)d2 (m/volt)L (m)
Value 1000 0.25 1.5 0.7 0.9 100 0.4 0.4 120 × 10−12 120 × 10−12 0.1 
For this specific WM, the open-loop compliance s 11 v and specific volume s 22 v are given by:
C o m p l i a n c e = s 11 v = ( 62.50 s 6 + 531.25 s 4 + 596.25 s 2 90 ) 112 , 500 s 8 + 468 , 250 s 6 + 551 , 550 s 4 + 153 , 450 s 2 + 3 , 024
and
S p e c i f i c V o l u m e = s 22 v = 125 s 2 ( 50 s 6 + 425 s 4 + 477 s 2 72 ) 112 , 500 s 8 + 468 , 250 s 6 + 551 , 550 s 4 + 153 , 450 s 2 + 3 , 024 .

Note that the poles of the open-loop compliance and specific volume are the same and are as given in Table II.

TABLE II.

Poles of the open-loop compliance and specific volume of the Willis material.

Poles1,23,45,67,8
Value ±1.5004i ±1.233i ±0.6069i ±0.1460i 
Poles1,23,45,67,8
Value ±1.5004i ±1.233i ±0.6069i ±0.1460i 

It is important to mention that all the poles are imaginary as expected as the AWM is undamped and the material, in its open-loop form, is marginally stable. This justifies considering the use of the proposed proportional-integral-derivative (PID) controllers described in Eqs. (37) and (47).

The optimization problem given by Eq. (40) is used to extract the optimal gain K υ with its components k p υ , k ¯ υ i , k ¯ υ d , a n d g r υ. To simplify the search for the optimal gain g r υ is taken = 1, i.e., the two piezo-actuators are equally controlled.

The search for the optimal gain is carried out, using a MATLAB Optimization Toolbox, with lower and upper bounds of the three design parameters k p υ , k ¯ υ i , a n d k ¯ υ d of 0.1 and 100. The search is initiated from an initial guess of (10,10,10). The optimal values attained by the search are k p υ = 0.84 , k ¯ υ i = 0.38 , a n d k ¯ υ d = 8.4 as displayed in Fig. 6(a). These values are reached after 20 iterations as shown in Fig. 6(b).

FIG. 6.

Optimal parameters of the controller k p υ , k ¯ υ i , a n d k ¯ υ d.

FIG. 6.

Optimal parameters of the controller k p υ , k ¯ υ i , a n d k ¯ υ d.

Close modal

The optimal gains result in closed-loop poles as listed in Table III.

TABLE III.

Poles of the closed-loop specific volume of the Willis material.

Poles12,34,56,78,910,11
Value ±0.6445i ±1.233i ±1.5493i −0.0002 ± 1.0013i −0.0594 ± 0.2029i 
Poles12,34,56,78,910,11
Value ±0.6445i ±1.233i ±1.5493i −0.0002 ± 1.0013i −0.0594 ± 0.2029i 

It can be seen that the first pole defines the integral component of the controller, poles 2 through 7 are complex-conjugate imaginary poles that do not destabilize the system, whereas poles 8 through 11 have damping components resulting from the damping effect generated by the PID controller.

With these inherent features, the closed-loop root locus plot shown in Fig. 7 emphasizes the fact that the optimal controller is stable as all the root loci lie in the left-hand-side of the Laplace s-plane.

FIG. 7.

Root-locus plot of the optimal closed-loop specific volume controller.

FIG. 7.

Root-locus plot of the optimal closed-loop specific volume controller.

Close modal

The behavior of the optimal specific volume controller is determined using the MATLAB-Simulink block diagram shown in Fig. 8.

FIG. 8.

MATLAB–Simulink block diagram of the specific volume controller.

FIG. 8.

MATLAB–Simulink block diagram of the specific volume controller.

Close modal

Figure 9 displays the time response of the controller when a specific volume of υR = 1 m3/kg is targeted with and without the effect of the external disturbance. Note that the effect of the external disturbance is exaggerated considerably by using the value of σ/p = 1 × 106 N s/(kg m3).

FIG. 9.

Time response of the specific volume controller for different levels of external disturbances.

FIG. 9.

Time response of the specific volume controller for different levels of external disturbances.

Close modal

The displayed results clearly demonstrate the effectiveness of the controller in targeting the desired physical property while rejecting the effect of the external disturbances.

Figure 10 displays the time response of the controller while targeting different desired specific volumes of υR from 0.1 to 10 m3/kg. The displayed results are obtained for a very high level of external disturbance of σ/p = 1 × 106 N s/(kg m3). It is evident that the controller has achieved the desired performance.

FIG. 10.

Time response of the specific volume controller for different desired reference specific volumes while rejecting the effect of external disturbances of σ/p = 1 × 106 N s/(kg m3).

FIG. 10.

Time response of the specific volume controller for different desired reference specific volumes while rejecting the effect of external disturbances of σ/p = 1 × 106 N s/(kg m3).

Close modal

The optimization problem given by Eq. (50) is used to extract the optimal gain K s with its components k p s , k ¯ s i , k ¯ s d , a n d g r s. To simplify the search for the optimal gain, g r s is taken equal to1. This means that the two piezo-actuators are equally controlled.

The search for the optimal gain is carried out, using a MATLAB Optimization Toolbox, with lower and upper bounds of the three design parameters k p s , k ¯ s i , a n d k ¯ s d of 0.01 and 100. The search is initiated from an initial guess of (50,50,1). The optimal values attained by the search are k p s = 55.84 , k ¯ s i = 0.01 , and k ¯ s d = 0.664 as displayed in Fig. 11(a). These values are reached after 72 iterations as shown in Fig. 11(b).

FIG. 11.

Optimal parameters of the compliance controller k p s , k ¯ s i , a n d k ¯ s d.

FIG. 11.

Optimal parameters of the compliance controller k p s , k ¯ s i , a n d k ¯ s d.

Close modal

The optimal gains result in closed-loop poles as listed in Table IV.

TABLE IV.

Poles of the closed-loop compliance of the Willis material.

Poles12,34,56,78,910,11
Value −1.668 × 105 ±1.235i −0.0004 ± 1.5493i −0.0056 ± 1.0000i −0.0022 ± 0.8207i −0.0001 ± 0.3111i 
Poles12,34,56,78,910,11
Value −1.668 × 105 ±1.235i −0.0004 ± 1.5493i −0.0056 ± 1.0000i −0.0022 ± 0.8207i −0.0001 ± 0.3111i 

It can be seen that the first pole defines the integral component of the controller, pole 1 is real negative, and poles 2 through 3 are complex-conjugate imaginary poles that do not destabilize the system, whereas poles 4 through 11 have damping components resulting from the damping effect generated by the PID controller.

With these inherent features, the closed-loop root locus plot shown in Fig. 12 emphasizes the fact that the optimal controller is stable as all the root loci lie in the left-hand-side of the Laplace s-plane.

FIG. 12.

Root-locus plot of the optimal closed-loop compliance controller.

FIG. 12.

Root-locus plot of the optimal closed-loop compliance controller.

Close modal

The behavior of the optimal compliance controller is determined using the MATLAB-Simulink block diagram shown in Fig. 13.

FIG. 13.

MATLAB–Simulink block diagram of the compliance controller.

FIG. 13.

MATLAB–Simulink block diagram of the compliance controller.

Close modal

Figure 14 displays the time response of the controller when a compliance of sR = 1 m2/N is targeted with and without the effect of the external disturbance. Note that the effect of the external disturbance is exaggerated considerably by using the value of p/σ = 1 × 106 kg m3/ s N.

FIG. 14.

Time response of the compliance controller for different levels of external disturbances.

FIG. 14.

Time response of the compliance controller for different levels of external disturbances.

Close modal

Figure 15 displays the time response of the compliance controller while targeting different desired compliances of sR from 0.1 to 10 m2/N. The displayed results are obtained for a very high level of external disturbance of p/σ = 1 × 106 kg m3/s N. It is evident that the controller has achieved the desired performance.

FIG. 15.

Time response of the compliance controller for different desired reference compliances while rejecting the effect of external disturbances of p/σ = 1 × 106 kg m3/s N.

FIG. 15.

Time response of the compliance controller for different desired reference compliances while rejecting the effect of external disturbances of p/σ = 1 × 106 kg m3/s N.

Close modal

This paper has presented a comprehensive modeling, analysis, design, and simulation of the characteristics of piezoelectric-based active Willis metamaterials (AWM). The developed model of the AWM are generated, from the first principles using the Lagrange dynamics formulation, to establish the constitutive coupling between the passive Willis metamaterial and active piezoelectric components. The generated form is compatible for direct application of the control systems theory. It therefore enables its use in actuation, sensing, and as an integrated controller that can be analyzed, designed, and optimized using the classical, optimal, and robust control system theories.

With this developed controlled-based structure of the AWM, the inherent and powerful capabilities of the AWM, which lie in its ability to robustly control the material properties themselves such as compliance (or stiffness) and specific volume (or density), are demonstrated in great detail via several numerical examples.

The presented examples utilize the developed theoretical models of the compliance (or stiffness) and specific volume (or density) to generate controlled-based structures that include reference command, system dynamics, external disturbances, and the controller. The gains of the controllers are designed to minimize the tracking errors and the disturbance errors and at the same time minimize the control effort using the H-2 minimization approach. In this manner, the tracking of the reference command (compliance or specific volume) is achieved while rejecting the effect of the external disturbance, accommodating any parameter uncertainty and ensuring a balance with the required control effort. Implementation of the controller is carried out using MATLAB-SIMULINK and is used to demonstrate the robustness of the controller in the presence of external disturbances and any parameter uncertainty.

Particular emphasis is placed on ensuring stability of the closed-loop system by computing the closed-loop poles and also using the root locus plot to visualize the stability in the Laplace domain.

It is important to note that although the presented study has focused on a one-dimensional modeling of the AWM, it can be easily extended to 2 and 3 dimensional arrays using similar components or, more generally, for other forms of continuum as plates and shells. Furthermore, accounting for the anisotropic properties of the piezoelectric ingredients can be easily included. All these additional features are a natural extension of the present work.

Work is also in progress to demonstrate the ability of the AWM is generating controllable bandgap and dispersion characteristics.

Furthermore, extensive effort is now in progress to demonstrate the controllable capabilities of the AWM experimentally and also to validate the developed models presented in this paper.

Finally, demonstrating the ability of the AWM in controlling the compliance (i.e. Bulk's modulus) and specific volume (i.e. density) provides various potential opportunities for using the AWM in numerous critical applications such as cloaking, beam shifting, and beam.

A

Cross sectional area of the piezoelectric material

Cs

Piezo-capacitance matrix

C p 1 , 2 s

Capacitances of the piezoelectric springs under strain-free conditions

d

Piezo-strain coefficients matrix

d

Piezoelectric electric strain coupling coefficients

d1,2

Piezoelectric electric strain coupling coefficients of the two springs

D

Piezoelectric electrical displacement

D

Parameter given by = [ μ 2 s ¯ 2 + ( 1 κ 1 2 ) + k r ( 1 κ 2 2 ) ]

E

Piezoelectric electric field

E1,2

Piezoelectric electric fields applied on the two springs

Es

Control effort associated with compliance

E s 2

H-2 norms of the control effort associated with compliance

Ev

Control effort associated with specific volume

E υ 2

H-2 norms of the control effort associated with specific volume

F

Force

FL,R Forces acting on the left and right masses of the unit cell of the Willis material

[G]

Control gain matrix

G 1 s

G 1 s = 1 / ( 1 + K s )

G 2 s

G 2 s = s 12 v / ( 1 + K s )

G 1 υ

G 1 υ = 1 / ( 1 + K υ )

G 2 υ

G 2 υ = s 21 v / ( 1 + K υ )

g

Charge sensitivity matrix

g1,2,v

Specific volume control gains

g3,4,s

Compliance volume control gains

g r s

g r s = g 4 g s

g r υ

g r υ = g 2 g υ

[H]

Matrix defined in Eq. (14)

H s 2

H-2 norms of compliance

H υ 2

H-2 norms of specific volume

h1,2

Piezoelectric coupling coefficients for the two springs

[KD]

Matrix defined in Eq. (14)

KE

Kinetic energy

k p υ

Proportional gain of specific volume

k p s

Proportional gain of compliance

Ks

K s = g s σ ¯ ( d 11 + g r s d 12 )

K υ = g υ p ¯ ( d 21 + g r υ d 22 )

k s d

= k p s k ¯ s d = D e r i v a t i v e g a i n

k s i

= k p s k ¯ s i = I n t e g r a l g a i n

k υ d

= k p υ k ¯ υ d = D e r i v a t i v e g a i n

k υ i

= k p υ k ¯ υ i = I n t e g r a l g a i n

k 1 , 2 D

Stiffness of the piezoelectric springs under open-circuit conditions

kr

The ratio k r = k 2 D / k 1 D

L

The distance between any of the two masses at equilibrium

L

Lagrangian

[M(1),(2)]

Stress–momentum transformation matrices defined in Eq. (16)

m1,2,3

Masses of the unit cell of the Willis material

PE

Potential energy

p

Momentum

p ¯

Nominal momentum

Q1,2

Electric charges generated by the two piezoelectric springs

S

Piezoelectric strain

s

Laplace complex number

sE

Short-circuit compliance of the piezoelectric material

sv

Short-circuit compliance and density matrix

s ¯

Dimensionless Laplace number = s / ω n

T

Piezoelectric stress

tp

Thickness of the piezoelectric material

uL,m,R

Deflections of the three masses of the unit cell of the Willis material

V1,2

Piezoelectric electric voltages applied on the two springs

v

Velocity

vR

Reference velocity

x

Deflection of the piezoelectric material

ɛ

Strain

ε R

Reference strain

ε s

Piezoelectric permittivity under strain-free conditions

ε T

Piezoelectric permittivity under stress-free conditions

κ 1 , 2 2

Electromechanical coupling factors of the two springs

μ2,3

The mass ratios m2/m1 and m3/m1

ω n 2

Natural frequency = k 1 D / m 1

σ

Stress

σ ¯

Nominal stress

υ

Specific volume

υ R

Reference specific volume

The author has no conflicts to disclose.

A. Baz: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The datasets generated and supporting the findings of this article are available from the corresponding author upon reasonable request.

A. Basic properties

The constitutive equations describe the interaction between the mechanical and electrical characteristics of piezoelectric films. In the one-dimensional case, these equations, as given in the IEEE Standard on Piezoelectricity,29 are
S = s E T + d E
(A1)
and
D = d T + ε T E .
(A2)

Equation (A1) describes the mechanical strain S due to a mechanical stress T and an electrical field E. In Eq. (A2), the electric displacement D is generated by the mechanical stress T and the electric field E. In these equations, sE, d, and ɛT denote the short-circuit compliance, the piezo-strain coefficient, and stress-free permittivity, respectively.

Note that the stress T, strain S, electrical displacement D, and electrical potential E can be expressed in terms of the actuator physical and geometrical parameters as follows:
T = F A , S = x t p , D = Q A , and E = V t p ,
(A3)
where F, A, x, tp, and Q denote the force, cross sectional area, displacement, thickness, and electrical charge of the piezo-spring.
Combining Eqs. (A3), (A1) and (A2) yields
F = k p D x h Q
(A4)
and
V = h x + 1 C p s Q ,
(A5)
where k p D = c D A t p denotes the open-circuit mechanical stiffness, c D = 1 s E ( 1 κ 2 ) , κ 2 = d 2 ε T s E is the “Electro-Mechanical Coupling Factor” (EMCF), and C p s = A ε S t p denotes the strain-free electrical capacitance of the piezo-spring. Also, ε s = ε T ( 1 κ 2 ) is the strain-free permittivity, and h = κ 2 d ( 1 κ 2 ) .

B. Potential energy

The potential energy (PE) can be determined from Eqs. (A4) and (A5) as follows:
PE = F d x + V d Q = ( k p D x ) d x + ( 1 C p s Q ) d Q h ( Q d x + x d Q ) = 1 2 k p D x 2 + 1 2 1 C p s Q 2 h Q x .
(A6)
1.
J. R.
Willis
, “
Variational principles for dynamic problems for inhomogeneous elastic media
,”
Wave Motion
3
(
1
),
1
11
(
1981
).
2.
M. B.
Muhlestein
, “Willis coupling in acoustic and elastic metamaterials,” Ph.D. thesis (The University of Texas at Austin, 2016).
3.
M. B.
Muhlestein
,
C. F.
Sieck
,
P. S.
Wilson
, and
M. R.
Haberman
, “
Experimental evidence of Willis coupling in a one-dimensional effective material element
,”
Nat. Commun.
8
,
15625
(
2017
).
4.
A.
Melnikov
,
Y. K.
Chiang
,
L.
Quan
,
S.
Oberst
,
A.
Alù
,
S.
Marburg
, and
D.
Powell
, “
Acoustic meta-atom with experimentally verified maximum Willis coupling
,”
Nat. Commun.
10
,
3148
(
2019
).
5.
J.-P.
Groby
,
M.
Malléjac
,
A.
Merkel
,
V.
Romero-García
,
V.
Tournat
,
D.
Torrent
, and
J.
Li
, “
Analytical modeling of one-dimensional resonant asymmetric and reciprocal acoustic structures as Willis materials
,”
New J. Phys.
23
(
5
),
053020
(
2021
).
6.
Y.
Liu
,
Z.
Liang
,
J.
Zhu
,
L.
Xia
,
O.
Mondain-Monval
,
T.
Brunet
,
A.
Alù
, and
J.
Li
, “
Willis metamaterial on a structured beam
,”
Phys. Rev. X
9
,
011040
(
2019
).
7.
G. W.
Milton
and
J. R.
Willis
, “
On modifications of Newton’s second law and linear continuum elastodynamics
,”
Proc. R. Soc. A
463
,
855
880
(
2007
).
8.
R.
Pernas-Salomón
and
G.
Shmuel
, “
Fundamental principles for generalized Willis metamaterials
,”
Phys. Rev. Appl.
14
,
064005
(
2020
).
9.
C. F.
Sieck
,
A.
Alù
, and
M. R.
Haberman
, “
Origins of Willis coupling and acoustic bianisotropy in acoustic metamaterials through source-driven homogenization
,”
Phys. Rev. B
96
,
104303
(
2017
).
10.
M. B.
Muhlestein
,
C. F.
Sieck
,
A.
Alù
, and
M. R.
Haberman
, “
Reciprocity, passivity and causality in Willis materials
,”
Proc. R. Soc. A
472
,
20160604
(
2016
).
11.
C.
Cho
,
X.
Wen
,
N.
Park
, and
J.
Li
, “
Acoustic Willis meta-atom beyond the bounds of passivity and reciprocity
,”
Commun. Phys.
4
,
82
(
2021
)..
12.
A.
Merkel
,
V.
Romero-García
,
J.-P.
Groby
,
J.
Li
, and
J.
Christensen
, “
Unidirectional zero sonic reflection in passive PT-symmetric Willis media
,”
Phys. Rev. B
98
,
201102
(
2018
).
13.
J.
Li
,
C.
Shen
,
A.
Díaz-Rubio
,
S.
Tretyakov
, and
S.
Cummer
, “
Systematic design and experimental demonstration of bianisotropic metasurfaces for scattering-free manipulation of acoustic wavefronts
,”
Nat. Commun.
9
,
1342
(
2018
).
14.
L.
Quan
,
D. L.
Sounas
, and
A.
Alù
, “
Nonreciprocal Willis coupling in zero-index moving media
,”
Phys. Rev. Lett.
123
,
064301
(
2019
).
15.
C.
Olivier
,
G.
Poignand
,
M.
Malléjac
,
V.
Romero-García
,
G.
Penelet
,
A.
Merkel
,
D.
Torrent
,
J.
Li
,
J.
Christensen
, and
J.-P.
Groby
, “
Nonreciprocal and even Willis couplings in periodic thermoacoustic amplifiers
,”
Phys. Rev. B
104
,
184109
(
2021
).
16.
R.
Pernas-Salomón
and
G.
Shmuel
, “
Symmetry breaking creates electro-momentum coupling in piezoelectric metamaterials
,”
J. Mech. Phys. Solids
134
,
103770
(
2020
).
17.
Y.
Zhai
,
H.-S.
Kwon
, and
B.-I.
Popa
, “
Active Willis metamaterials for ultracompact nonreciprocal linear acoustic devices
,”
Phys. Rev. B
99
,
220301
(
2019
).
18.
K.
Muhafra
,
M. R.
Haberman
, and
G.
Shmuel
, “
Discrete one-dimensional models for the electromomentum coupling
,”
Phys. Rev. Appl.
20
,
014042
(
2023
).
19.
Y.
Zhai
, “An approach for active acoustic metamaterial design implementation and analysis,” Ph.D. thesis (University of Michigan, 2022).
20.
S. R.
Craig
,
B.
Wang
,
X.
Su
,
D.
Banerjee
,
P. J.
Welch
,
M. C.
Yip
,
Y.
Hu
, and
C.
Shi
, “
Extreme material parameters accessible by active acoustic metamaterials with Willis coupling
,”
J. Acoust. Soc. Am.
151
,
1722
1729
(
2022
).
21.
M.
Kosta
,
A.
Muhafra
,
R.
Pernas-Salómon
,
G.
Shmuel
, and
O.
Amir
, “
Maximizing the electro-momentum coupling in piezoelectric laminates
,”
Int. J. Solids Struct.
254–255
,
111909
(
2022
).
22.
A.
Muhafra
,
M.
Kosta
,
D.
Torrent
,
R.
Pernas-Salomón
, and
G.
Shmuel
, “
Homogenization of piezoelectric planar Willis materials undergoing antiplane shear
,”
Wave Motion
108
,
102833
(
2022
).
23.
Z.
Zhang
,
J.-H.
Lee
, and
G. X.
Gu
, “
Rational design of piezoelectric metamaterials with tailored electro-momentum coupling
,”
Extreme Mech. Lett.
55
,
101785
(
2022
).
24.
H. D.
Huynh
,
X.
Zhuang
,
H. S.
Park
,
S. S.
Nanthakumar
,
Y.
Jin
, and
T.
Rabczuk
, “
Maximizing electro-momentum coupling in generalized 2D Willis metamaterials
,”
Extreme Mech. Lett.
61
,
101981
(
2023
).
25.
J.-H.
Lee
,
Z.
Zhang
, and
G. X.
Gu
, “
Reaching new levels of wave scattering via piezoelectric metamaterials and electro-momentum coupling
,”
J. Acoust. Soc. Am.
153
,
A163
(
2023
).
26.
A.
Baz
, “
Active acoustic metamaterial with tunable effective density using a disturbance rejection controller
,”
J. Appl. Phys.
125
(
7
),
074503
(
2019
).
27.
A.
Baz
, “Active acoustic metamaterials with programmable densities using an H-∞ controller,” in
Proceedings of the ASME 2018 International Mechanical Engineering Congress and Exposition
, V011T01A014, IMECE2018-87749 (ASME, 2019).
28.
A.
Baz
,
Smart Structures Concepts and Applications
(
World Scientific
,
2024
).
29.
ANSI/IEEE STD 176-1987, IEEE Standard on Piezoelectricity, The Institute of Electrical and Electronics Engineers, Inc., 1987.