## Abstract

The effect of microstructural properties on the wave dispersion in linear elastic membranes is addressed in this paper. A periodic spring-mass lattice at the lower level of observation is continualized and a gradient-enriched membrane model is obtained to account for the characteristic microstructural length scale of the material. In the first part of this study, analytical investigations show that the proposed model is able to correctly capture the physical phenomena of wave dispersion in microstructured membrane which is overlooked by classical continuum theories. In the second part, a finite-element discretization of microstructured membrane is formulated by introducing the pertinent inertia and stiffness terms. Importantly, the proposed modifications do not increase the size of the problem compared wiith classical elasticity. Numerical simulations confirm that the vibrational properties are affected by the microstructural characteristics of the material, particularly in the high-frequency regime.

## 1. Introduction

A tensioned membrane can be considered as the simplest example of a two-dimensional vibrating system, although more complex structures may also be considered as membranes, the fundamental characteristic of which is the incapability of conveying moments or shear forces. Currently, innovative practical applications exist for these structural models, such as antennas and solar sails for space applications, biotechnological prostheses and lightweight roofs designed in architecture and structural engineering (Jenkins & Korde 2006). The use of composite materials offers excellent mechanical performance combined with lightweight characteristics appropriate for membrane structures. Composite membranes tend to possess relevant microstructure that often significantly influences the corresponding macroscopic response. Consequently, accurate modelling of such structures and materials requires that the microstructure is accounted for.

In mechanical problems involving materials with microstructure, the existence of microscopic characteristic lengths, whose origin may be the lattice period of a discrete system or the inclusion size of an inhomogeneous material, introduces a size dependence into the structural response on the macroscopic level. These effects may be accurately analysed, including each individual microstructural element in a discrete model. However, this approach leads to a considerable computational effort because the number of degrees of freedom is normally very large. The strategy of replacing such a discrete model with an equivalent continuum, characterizing the average mechanical behaviour of the actual discrete or inhomogeneous material, offers an efficient alternative, as the continuum mechanics approach is computationally less expensive.

Although classical continuum theories are adequate for many static loading conditions, they may not be appropriate for problems involving dynamic analysis in the high-frequency regime, wherein microstructural effects are more significant. For instance, classical continuum theories cannot capture wave dispersion. Dispersion is the phenomenon that waves with different wavelengths propagate with different velocities, and wave dispersion in heterogeneous or discrete media is extensively documented by experimentalists (Erofeyev 2003; Jakata & Every 2008).

Waves that propagate through a classical medium are not dispersive, i.e. the wave speed is constant, independent of wavelength. In contrast, enhanced continuum models account for the intrinsic length of microstructural details in the macrostructural material model, normally by means of additional (temporal or spatial) derivatives of the relevant state variables. It was shown that the wave dispersion can be captured by the addition of higher-order inertia terms and higher-order stiffness terms in the governing equations, as discussed in the recent overview by Papargyri-Beskou *et al.* (2009). To obtain a transparent link between the microstructural material properties and the higher-order terms on the macroscopic level, continualization or homogenization strategies were proposed to account for the heterogeneous or discrete nature of materials—in the context of wave propagation studies, see, for instance, the studies of Chang & Gao (1995), Mühlhaus & Oka (1996), Chen & Fish (2001), Suiker *et al.* (2001), Fish *et al.* (2002), Metrikine & Askes (2002) or Suiker & de Borst (2005).

Finite-element implementations are essential for widespread dissemination of material models. Unfortunately, finite-element implementations tend to be complicated for the enhanced continuum models mentioned above. As discussed, for example, by Askes *et al.* (2008), many formulations of enhanced continua require special treatment, such as Hermitian or mixed formulations, which makes their implementation in available finite-element packages cumbersome. Thus, it is of interest to pursue enhanced continuum formulations that allow straightforward finite-element implementation.

This contribution addresses a number of aspects of enhanced continuum models mentioned above. The main aim is to validate the accuracy of enhanced continuum analysis of a membrane with microstructure subject to dynamic loading. The continualization approach of Rosenau (1987) and Andrianov & Awrejcewicz (2008) will be applied in §2 to translate microstructural properties into a macroscopic model. An analysis of the dispersion properties of transverse waves propagating along different directions is presented in §3. The spatial discretization with standard finite elements is detailed in §4. Finally, numerical analyses of a microstructured membrane are treated in §5 with the aim to verify that the enhanced continuum model provides more realistic results than the classical continuum theory.

## 2. Transverse vibrations of a two-dimensional lattice and continualization methods

In this section, the equations governing the transverse vibrations of a linear elastic two-dimensional lattice have been derived. The continualization approach is employed for the construction of continuum models from associated periodic structures. Truncated Taylor series expansions are used to approximate the difference equations by differential equations. In this manner, the microstructural effects are accounted for by means of higher-order gradients naturally appearing in the continuum equations. The equation of motion obtained via the standard continualization approach is shown to reduce to classical elastic continuum when only the first terms in the series are accounted for. Higher-order terms increase the accuracy of the approximation, but also the order of the partial differential equations, which requires supplementary boundary conditions to define a well-posed problem. In order to avoid this inconvenience, the Padé approximants are introduced, yielding a more effective continuum model; incidentally, Padé approximants also help in overcoming instabilities as explained in §3.

### (a) Discrete model

The two-dimensional spring-mass lattice that will be analysed consists of a repeating arrangement of *N*_{1}×*N*_{2} identical particles of mass *M*, connected by pretensioned springs. The pretensioning of the springs occurs entirely in plane, but provides stiffness against out-of-plane deflection. The interparticle distances in the *x* and *y* directions are denoted by ℓ_{x} and ℓ_{y}, respectively. As outlined in Rosenau (1987), owing to the spatial periodicity of the system, the dynamics of such structures can be studied by considering only a unit cell, as depicted in figure 1. If the lattice experiences small displacements in the transverse (i.e. *z*) direction compared with its initial position in the *x*–*y* plane, the membrane vibrations can be studied in the linear regime (Andrianov & Awrejcewicz 2008).

The equation of motion governing the free transverse vibrations of an element located in the *m*th row and *n*th column is obtained by balance of forces in transversal direction as
2.1
where *u*_{m,n}=*u*_{m,n}(*t*)=*u*(*x*_{m},*y*_{n};*t*) is the transverse displacement of the mass and *F*_{i,j} is the recall force of the stretched springs connecting the nearest-neighbour particles.

For small displacements, we may assume that the contribution of *F*_{i,j} at all points in the lattice remains close to its equilibrium value, and each contribution has the form
2.2a
and
2.2b
where now *F*_{i} is the force of the spring in direction *i*.

For the sake of simplicity, we employ a square lattice, so that *F*_{x}=*F*_{y} and ℓ_{x}=ℓ_{y}, and, inserting equations (2.2) into equation (2.1), the equation of motion reads
2.3
where *T*=(*F*/ℓ) is an in-plane tensile load. The following boundary and initial conditions are applied:
2.4a
and
2.4b
To continualize equation (2.3), the simplest strategy suggests to replace the discrete degrees of freedom by a continuous variable field.

### (b) Classical continuum approximation

Following the hypothesis of a dense lattice, a continuum approximation involves introducing a continuum displacement field *u*_{m,n}≡*u*(*x*,*y*;*t*), whereas the displacements of the neighbouring particles *u*_{m±1,n} and *u*_{m,n±1} are replaced by *u*(*x*±ℓ,*y*) and *u*(*x*,*y*±ℓ), respectively. Using Taylor series expansions, the displacements of the neighbouring particles can be approximated according to
2.5a
and
2.5b

For notational simplicity, the dependence of *u* on the spatial coordinates and on the time *t* is not explicitly shown. The accuracy of the derived continuum depends on the number of terms included in the Taylor series. Taking into account the terms up to the second derivatives and introducing in equation (2.3) the continuum field variables, a membrane-like continuum approximation is obtained,
2.6
where *ρ*=(*M*/ℓ^{2}) is the mass density. If the uniform tension *T* is sufficiently large, the transverse displacement is small (so that *T* does not change significantly). Equation (2.6) represents the classical equation of transverse motion for a vibrating membrane uniformly stretched, bounded by a surface *Γ*, with initial conditions *u*_{0} and (∂*u*_{0}/∂*t*), and boundary conditions fixed along a portion of the edge *Γ*_{D}, i.e. *u*=0 on *Γ*_{D} and free along the boundary *Γ*_{N}, i.e. *T*(∂*u*/∂*n*)=0 on *Γ*_{N}. The applicability of classical continuum models in dynamics is, however, restricted to a limited range of frequencies as will be highlighted in §3.

It is also worth noting that, for the sake of simplicity, coupling between transversal and longitudinal vibration of the membrane is neglected in the proposed formulation. This simplification is acceptable as long as the slope of membrane is small. Further studies will be needed to assess the effect of the coupled motion of the membrane.

### (c) The quasi-continuum method and Padé approximation

Although the application of a model with derivative of more than second order allows to include the effect of microstructure via the parameter ℓ in the continuum model, the higher order of the resulting differential equation complicates the solution of continuous problems. An alternative continualization procedure can be realized using the quasi-continuum approximations. This technique consists of the reduction of the differential-difference equation (2.3) to a partial differential equation. The continuum model is not obtained by simply replacing the difference operator in equation (2.3) by the analogous differential operator. In the first stage, the backward and forward shift operators are introduced, and for the displacements of adjacent masses, the following symbolic notations hold:
2.7a
and
2.7b
By means of equations (2.7*a*,*b*) and the relation *u*_{m,n}=*u*, the discrete equation (2.3) is converted into a pseudo-differential equation as follows:
2.8
The differential operators in square brackets can be developed into Taylor series in the neighbourhood of zero (McLaurin series),
2.9a
and
2.9b
Taking into account only the first term in the expansions (2.9), one obtains again the classical continuum approximation given in equation (2.6). Furthermore, the accuracy of the approximation can be increased keeping the next terms in the series, and the higher-order approximation reads
2.10
where ∇^{2}=(∂^{2}/∂*x*^{2})+(∂^{2}/∂*y*^{2}).

An alternative continuum model has been developed by using the one-point Padé approximation of the differential operator appearing on the right-hand side of equation (2.10). For more details regarding the Padé approximants, we refer the reader to the work of Rosenau (1986), Wattis (2000), Kevrekidis *et al.* (2002) and Andrianov & Awrejcewicz (2005). Taking into account only two terms in the expansion, the operator in square brackets can be developed by means of the Padé approximation
2.11
which leads to the following continuum model:
2.12
The characteristic length ℓ accompanies the higher-order inertia term on the left-hand side as well as the higher-order stiffness term on the right-hand side.

## 3. Dispersion analysis

The dispersive character of plane waves owing to inherent material characteristic length is consistent with experimental observations. Specifically, this effect becomes significant when the wavelength of the physical phenomena of interest decreases to the order of the internal length, as, for example, in problems related to shock waves. The regularization procedures introduced in the previous section are developed with the aim to define a non-local continuum that models dispersion phenomena by capturing the characteristics of the microstructure. The effectiveness of the enhanced continuum models to represent the wave propagation characteristics of a lattice system can be assessed by evaluating the dispersion relations.

### (a) Dispersion relations for the lattice model and the classical continuum approximation

In the analysis of plane-wave propagation in two-dimensional periodic discrete structures, a dispersion equation can be obtained by substituting a harmonic wave into equation (2.3) as
3.1
where *ω* is the angular frequency, *A* is the amplitude, is the imaginary unit and *x*_{m}=*m*ℓ and *y*_{n}=*n*ℓ are the spatial coordinates, respectively. In a two-dimensional model, transverse waves travel along an arbitrary direction that is at an angle *θ*, which represents the angle between the wavevector ** k** and the

*x*axis. In equation (3.1), and define the components of the wavenumber in the

*x*and

*y*directions, respectively. The dispersion relation is a function

*ω*=

*f*(

*k*

_{x},

*k*

_{y}), which contains information of the wave propagation characteristics of the considered material and, from the equation of motion (2.3), for the discrete lattice model is obtained as follows: 3.2 where is a non-dimensional frequency parameter. In the analysis of wave dispersion, the non-dimensional wavenumbers

*k*

_{x}ℓ and

*k*

_{y}ℓ are varied to investigate the frequencies at which waves propagate and the direction of the body wave propagation.

The dispersion relation expressed in equation (3.2) is evidently periodic in the interval *k*_{x}ℓ,*k*_{y}ℓ∈[−*π*,*π*], and it is also symmetric with respect to the axes *k*_{x}ℓ=0 and *k*_{y}ℓ=0. These properties allow us to restrict the analysis to the range *k*_{x}ℓ,*k*_{y}ℓ∈[0,*π*] that reflects the so-called first Brillouin zone (Brillouin 1946). A particularly convenient representation of the dispersion relations employs iso-frequency curves whose graphical interpretation provides information on the dispersion properties. The curves determined by equation (3.2) are represented in normalized form (*k*_{x}ℓ−*k*_{y}ℓ plane for fixed values of ) in figure 2*a*.

The wave group velocity *c*_{g} for undamped structures is the velocity of propagation of vibrational energy, and may be expressed in the form
3.3
where ** i** and

**are the unit vectors in the**

*j**x*and

*y*directions, respectively. It should be noted that the group velocity is normal to the curves in the

*k*

_{x}−

*k*

_{y}plane for fixed

*ω*and its magnitude is equal to the gradient of the dispersion surface. Hence, the direction of propagation of a wave at a given frequency can be estimated by taking the normal to the iso-frequency contour lines (Langley 1997). This property, in particular, can be used to identify regions within the structure where waves at certain frequencies do not propagate in specified directions.

In §2, the discrete equation of motion has been transformed in a continuum model taking advantage of the continualization approach. As previously illustrated, if, in the series expansion, only the terms up to the second derivatives are retained, the equation of motion yields the wave equation based upon a classical continuum model expressed by equation (2.6). Substituting into the equation of motion the harmonic wave of the continuum form
3.4
we find the solution
3.5
This equation shows that the classical continuum is not dispersive because the frequency *ω* is directly proportional to the wavenumber *k*.

The dispersion relations of the square lattice and the classical continuum model, as represented by equations (3.2) and (3.5), respectively, have been plotted in figure 2*b*, which displays the behaviour of the dimensionless frequency in terms of the normalized wavenumber *k*ℓ. The dispersion curve for the continuum model is represented as a solid straight line through the origin, while the discontinuous lines represent the dispersion properties of the lattice for different angles *θ* in the range [0,*π*/4]. The comparison illustrates that the physical motion of the lattice system is well described by the classical continuum model for longer wave limits. This means that when the wavelengths are much larger than the structural inhomogeneity of the lattice ( or *k*ℓ≪1), the waves are not affected by the heterogeneity of the membrane.

On the other hand, the dispersive effects become evident in the short-wave range, when the phase velocity *c*_{p} of the harmonic waves, that is, the speed at which the phase of any one component of the wave travels, is different from the group velocity *c*_{g}. The phase velocity *c*_{p} can be expressed as
3.6
and is graphically represented by the secant slope of the dispersion curves, while the group velocity *c*_{g} is represented by the tangential slope. By inspection of the curves, the local continuum model appears unable to predict the decrease of phase velocity for increasing wavenumbers.

### (b) Dispersion relations for higher-order continua

Retaining more terms in the series gives an enhanced continuum model represented by equation (2.10). Explicitly, we seek solutions of the harmonic form expressed by equation (3.4) that, inserted into equation (2.10), gives the relationship between the frequency *ω* and the wavenumber *k*, which in non-dimensional form turns out to be
3.7
From a mathematical point of view, the model governed by equation (3.7) produces wave dispersion if the function under the square root is greater than zero. When *k* exceeds a critical value for a fixed direction of propagation, the wave does not possess real frequency and the solution is unstable. The dispersion curves for the higher-order continuum, depicted in figure 3, show a downward branch crossing the axis at values of *k*ℓ at which instability occurs. Because of the limitation owing to the stability criterion, the applicability of equation (2.10) is restricted to the dynamic phenomena that do not involve short waves.

If asymptotic analysis is resorted to, a way to deal with the instability problem in the fourth-order continuum model suggests replacing the fourth-order derivative in equation (2.10) with mixed double time double space derivatives (Pichugin *et al.* 2008). With the aim of restoring, in the continuum approximation, the dispersion property owing to the discreteness of the square lattice while, at the same, avoiding instabilities, the Padé approximations can be employed. In the following, the model introduced in equation (2.12) is referred to as a Padé continuum.

Solving the dispersion equation for the expected harmonic travelling wave, the following expression can be found: 3.8 It is worth remarking that the model is unconditionally stable, as the frequencies defined by equation (3.8) are real in the complete range of wavenumbers.

In figure 4, the dispersion curves for the Padé model are plotted in accordance with the relations obtained for the square lattice and the classical continuum model for different directions of propagation. The plots show that, at the long-range limit , the tangential slopes of the dispersion curves, corresponding to the group velocity *v*_{g}=(∂*ω*/∂*k*), are equal to those of the lattice model and the local continuum. For increasing wavenumbers, curves start to deviate from those of the discrete lattice model. However, they are located below the bold line that represents the local continuum model, and the group velocity is smaller than the one obtained for the local model, for all directions.

We note that the dispersion curve for the Padé continuum with *θ*=0 tends to a horizontal asymptote. In this case, the slope (group velocity) for transverse waves is zero and frequencies above the limit value cannot propagate. Thus, the model acts as a filter, where only relatively low frequencies are transmitted along the *θ*=0 direction.

By considering its definition, the group velocity of waves for the Padé continuum can be found to depend on the wavenumber and the direction of propagation as 3.9 where is the classical wave speed. This expression allows us to verify that propagating waves cannot transfer energy infinitely fast. In the long-wave limit, the group velocity takes the following form: 3.10 and for short waves 3.11 Thus, the model governed by equation (2.12) is stable and provides lower and upper bounds for the speed of energy transfer by propagating waves (Metrikine 2006).

## 4. Discretization of the equation of motion

In this section, a spatial discretization based on the finite-element method is carried out and the conditions for the convergence are reported.

### (a) Finite-element equations

A weak form of the initial-boundary value problem is taken by premultiplying equation (2.12) by a weight function *v* and integrating over the domain *Ω*, that is
4.1
The term *p* represents the influence of a transverse applied force per unit area, and it has been introduced in the last integral of equation (4.1) in order to generalize the initial boundary-value problem. Next, the divergence theorem is used to distribute the differentiation among *v* and *u*, so the weak form of equation (2.12) is written as
4.2
where *n*_{x} and *n*_{y} are the components of the unit normal vector on the boundary Γ, and a superimposed dot denotes the full derivative with respect to time. On the right-hand side of equation (4.2), the boundary integrals are collected whose contribution will be discussed afterwards in relation to the precise form of the boundary conditions of the problem. At this stage, we assume homogeneous boundary conditions, so the boundary integrals cancel.

In deriving the finite-element equations, the computational domain *Ω* is partitioned into bilinear isoparametric quadrilateral finite elements. This finite element is designated in the sequel as Q4 finite element (i.e. a quadrilateral finite element with four nodes). The classical nodal interpolation procedure of the displacement field is assumed, and the time dependence is separated from the spatial variation. The unknown displacement field *u*(*x*,*y*,*t*) is approximated over a typical finite element *Ω*^{e} by using standard polynomial shape functions *N*_{j} (*j*=1,…,4), as follows:
4.3
where *d*_{j} (*j*=1,…,4) are the nodal displacements of the *j*th node of the element. The weight function *v* is also discretized with the same shape functions. Because the weak form involves first-order and second-order mixed partial derivatives, continuous shape functions can be used, avoiding the cumbersome continuity in the numerical implementation.

Substituting the finite-element approximation in the weak form, the corresponding matrix form of equation (4.2) is given by 4.4 where 4.5 are the classical elasticity contributions to the element mass and stiffness matrices, respectively, while 4.6a and 4.6b are the higher-order contributions to the element mass and stiffness matrices, and 4.7 is the nodal force vector given by the externally applied pressure.

The equation of motion can be easily discretized and solved in the time domain by applying Newmark’s time-integration method. For an unconditionally stable solution, the constant average acceleration scheme is used.

### (b) Convergence and accuracy of the finite-element solution

To prove the convergence of the numerical solution to the exact solution of the governing equations as the number of elements increases, a measurement for its quality is required. If an exact solution is not available, the convergence can be measured only by the fact that some conditions contained in the mathematical model must be satisfied at convergence (Reddy 1993).

The accuracy of the finite-element method for the membrane problem when the enhanced continuum model is used may be assessed, for instance, from the convergence of the strain energy. The strain energy actually contained in the finite-element model can be evaluated via
4.8
In order to illustrate the performance of the Q4 membrane elements, the convergence in strain energy is shown in figure 5*a*,*b* for two sequences of uniform and skewed meshes obtained by halving the element size for the static problem of a membrane simply supported on the boundary and loaded by a uniform pressure. The plots give an error estimate if the convergence of strain energy has been achieved. As both meshes are refined, the strain energy approaches the exact solution from below.

Another, more frequently used, approach of quantifying the discretization error is the relative energy norm defined as follows:
4.9
where the reference energy *U*_{ref} corresponds to the finest mesh. This measure of error is closely related to the error on the derivative of the computed displacement field.

Error estimates of the type (equation (4.9)) are very useful because they give qualitative information on the accuracy of the approximate solution, whether or not the true solution is known. The energy norm as a function of mesh refinement is reported in figure 5*c*,*d* in logarithmic scaling. The error in the energy norm is proportional to *h*, which represents a typical length of an element side, and is expressed as
4.10
where *k* is the degree at which the interpolant polynomial is complete and *m* is the order of the highest derivative of the displacement field in the weak form (Zienkiewicz & Taylor 2000; Akin 2005). The element investigated here uses a polynomial displacement complete to the degree *k*=1, and in the weak form, the highest derivative order is of first order in the spatial dimensions (*m*=1). Thus, a rate of convergence *p*=*k*+1−*m*=1 is expected. As can be observed from figure 5*c*,*d*, the convergence obtained in the numerical simulation compares very well with the theoretical value. Moreover, having compared the numerical results obtained for the two test meshes, it is evident that the convergence rate is not affected by the element distortion.

## 5. Numerical examples

In §2, the equation of motion of a microstructured membrane was derived by means of non-local continuum mechanics models, and in §4, the discretization of the equation in the spatial and temporal domain was presented. Based on the formulations obtained above, to examine the microstructure effects on the dynamic response of a membrane simply supported at the edges, the results obtained from classical and higher-order elasticity (Padé continuum) are discussed in this section.

The geometry and the boundary conditions of the squared microstructured membrane are given in figure 6*a*. The length and the width of the panel is 2*L*=2 m, the material density is *ρ*=1 kg m^{−2}. The parameter ℓ is the microstructural characteristic length. The restoring force in the vibrating membrane arises from the tension *T* at which it is stretched and the in-plane tension *T*=16 N m^{−1} is chosen for the current analysis.

Because the dynamical response of a structure is mainly dependent on the dynamic load applied and on the dynamic characteristic of the structure itself, the analysis of natural frequencies is an important topic. The modal properties of a membrane can be obtained by solving the standard matrix eigenproblem
5.1
where *ω* is the finite-element predicted natural frequency of the system and **d** is the corresponding natural mode of the system.

As analytical solutions are only available for the classical continuum model and a mesh refinement allows to obtain more accurate results, a convergence analysis has been carried out in order to assess how many Q4 elements should be used in the numerical calculations. The frequencies based on classic elasticity and higher-order models are given in figure 7. The results associated with ℓ=0 correspond to those of the classical elasticity theory where the microstructural effect is ignored. It can be seen that the difference of the frequencies for different values of ℓ becomes evident at higher vibrational modes, although this difference is negligible at lower vibrational mode numbers. Further, the frequency decreases when the characteristic internal length ℓ increases. Hence, the frequencies are always overestimated by the classical continuum model.

The second numerical application deals with the transient response of the microstructured membrane subjected at its centre to an impact excitation. The dynamic load has been modelled as an initial non-zero velocity *v*(*t*_{0})=1 m s^{−1} at the impact point. Owing to the biaxial symmetry of the problem, we use a uniform mesh of Q4 elements to discretize only the top-right quadrant of the domain. The axes of symmetry *x*=0 and *y*=0 become a portion of the boundary *Γ* of the computational domain and the kinematic considerations reported in figure 6 allow us to justify the assumption of homogeneous boundary conditions introduced in §4. The Newmark scheme is used for the time integration of the equations using the recommendations on time-step size proposed by Bennett & Askes (2009).

Figure 8 shows the wavefronts at times *t*=0.10, 0.20 and 0.25 s for classical elasticity, obtained taking ℓ=0, and the Padé model for microstructural length ℓ=0.01, respectively. As previously predicted in §3 by means of the dispersion analysis, the elastic waves obtained using the local elasticity theory propagate faster than the case of the enhanced continua. The wavefronts for both models are governed by the components with higher wavelengths (small *k*) because they travel faster than the shorter ones. In fact, the velocity of the long waves is not influenced by the direction of propagation, and the fronts have quite a circular shape. The directional properties of the microstructured membrane expected for the enhanced continuum model can be observed from figure 8*b*, focussing the attention on the short waves that propagate slower and whose velocities vary with the direction of propagation. The directional property of the Padé model can be better observed in figure 9, where elastic waves propagating from the excitation point to the boundary for different characteristic lengths and along two specific directions are plotted at time *t*=0.20 s.

## 6. Conclusions

This contribution addresses the simulation of wave dispersion in composite membranes with microstructure. As an alternative to the direct analysis of the discrete material that, at microscopic level, has been assumed here as a periodic elastic lattice, a continuum model approximating the actual material behaviour has been employed. A higher-order theory has been successfully resorted to, overcoming the limits of classical elasticity in the high-frequency regime.

In the first part of this study, continuum models linked with the underlying material have been introduced, taking into account the discreteness of the microstructure. From the application of different continualization approaches, it has been shown that the higher-order theories improve the continuum description by including the microstructural effects. Specifically, enhanced continuum theories based on the use of Padé approximants introduce higher-order inertia and stiffness terms in classical models and improve the continuum description by including the length parameter characterizing the heterogeneity. The analysis of closed-form dispersion relations demonstrated that the enhanced continuum models deliver physically consistent results also in the limit of short waves and do not generate numerical instability.

Further, the field equations have been discretized in space by means of simple four-noded finite elements. In the implementation, continuity is completely avoided. In fact, the weak form of the governing equations involves first-order and second-order mixed derivatives of the unknown field, thus continuous shape functions are sufficient. Importantly, the higher-order terms do *not* lead to an increase in system size. The evaluation of the convergence rate confirms the appropriateness of the adopted finite-element discretization.

Finally, the analysis of a microstructured membrane has been conducted. The vibration frequencies have been estimated by solving the standard matrix eigenproblem for the classical and enhanced continuum models. It has been found that the frequencies are always overestimated by the classical continuum model. In order to illustrate the dispersion of the elastic waves, the response of the membrane to impulsive point load has been simulated. Numerical results are in agreement with the analytical dispersion analysis presented in this paper.

## Acknowledgements

We gratefully acknowledge financial support of the Engineering and Physical Sciences Research Council, contract number EP/D041368/1.

## Footnotes

- Received October 4, 2009.
- Accepted December 4, 2009.

- © 2010 The Royal Society