Microsoft Word - numero_29_art_31 G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 351 Focussed on: Computational Mechanics and Mechanics of Materials in Italy Mixed methods for viscoelastodynamics and topology optimization Giacomo Maurelli, Nadia Maini, Paolo Venini Department of Civil Engineering and Architecture, University of Pavia, Via Ferrata 3, 27100 Pavia, Italy paolo.venini@unipv.it ABSTRACT. A truly-mixed approach for the analysis of viscoelastic structures and continua is presented. An additive decomposition of the stress state into a viscoelastic part and a purely elastic one is introduced along with an Hellinger-Reissner variational principle wherein the stress represents the main variable of the formulation whereas the kinematic descriptor (that in the case at hand is the velocity field) acts as Lagrange multiplier. The resulting problem is a Differential Algebraic Equation (DAE) because of the need to introduce static Lagrange multipliers to comply with the Cauchy boundary condition on the stress. The associated eigenvalue problem is known in the literature as constrained eigenvalue problem and poses several difficulties for its solution that are addressed in the paper. The second part of the paper proposes a topology optimization approach for the rationale design of viscoelastic structures and continua. Details concerning density interpolation, compliance problems and eigenvalue-based objectives are given. Worked numerical examples are presented concerning both the dynamic analysis of viscoelastic structures and their topology optimization. KEYWORDS. Viscoelasticity; Mixed finite elements; Topology optimization. INTRODUCTION iscoelasticity is a constitutive feature of materials that finds applications in several areas such as structural engineering (concrete), biomedics (soft tissues) and also synthetic materials and polymers fabrication [1]. From a modeling viewpoint, the most classic approach to viscoelasticity is based on hereditary integral formulations that rely on the interpretation of viscoelastic materials as materials with memory: the current stress is based on a time integral of the strain history. Given such a complex constitutive response, the design of viscoelastic structures is quite a challenging task that calls for a robust numerical model and a rationale design approach. Within this scenario, object of this paper is the proposal of a truly-mixed finite element method for the analysis of viscoelastic structures coupled to a topology optimization approach with the ultimate goal of the development of an integrated analysis-design tool for viscoelastic systems. For this paper's sake, eigenvalue-based objectives shall be pursued within a few numerical applications that are based on the purely elastic cases discussed in [2]. As to the mixed finite-element formulation, the truly-mixed setting in a general Hellinger-Reissner framework has been selected mainly because of the need to pursue an accurate approximation of the stress field, feature that is not shared by any standard displacement-based approach. Furthermore, such a formulation is well known to pass the inf-sup condition [3] even in the limiting case of incompressible materials with no further tricks or enrichments and this could represent a V G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 352 key issue when modeling and designing rubber devices such as base isolators. The adopted formulation resembles the one presented in [4] in that an additive decomposition of the stress field is considered that finds its motivation if the adoption of a Maxwell-type phenomenological model. However, a weak enforcement of the stress tensor symmetry is considered in [4] as opposed to a strong approach adopted herein. The Arnold-Winther finite element is used as to the discretization of the stress field [5]. To the author's knowledge, it has never been used for the analysis of viscoelastic structures whereas a few applications in a purely elastic framework are available in the literature [6]. Optimal design with respect to eigenvalues is a topic that has received much attention in the last decades and [7, 8] and [9] may be cited among the most interesting contributions in this respect. However, optimal design of nonclassically damped viscoelastic structures is far from being a mature topic and further research seems to be in order. To gain insight into such a complex problem, a viscoelastic thin-beam truly-mixed formulation is introduced that has the merit to allow a deep understanding of the spectral properties of the discretized structure so as to end up with convincing optimal topologies. The first three eigenvalues shall be object of optimization and a strict relation between the optimal density distribution and the relevant eigenmode shapes clearly determined. Though preliminary, not many such results are available in the literature and should open the way to more complex results concerning two-dimensional systems. FORMULATION AND DISCRETIZATION (CONTINUUM CASE) oal of this section is the definition of a truly-mixed variational formulation for two-dimensional viscoelastic continua. Based on a parallel phenomenological model, compatibility, equilibrium and viscoelastic laws are written in strong form so as to allow the definition of a truly mixed variational formulation of Hellinger-Reissner type. Eventually, the Arnold-Winther finite-element is introduced along with some technical details concerning its implementation. Figure 1: Standard solid phenomenological model. Strong form Reference is made to Fig. 1 for the standard viscoelastic solid model that is the basis for the continuum and thin-beam models to be developed hereafter. As usual when adopting Hellinger-Reissner variational principles, compliance tensors relating strains to stresses are introduced that allow one to write 0 0 0 0 1 1 ( ) ( ) E V E A A v A v               (1) where 0EA and 0 VA are the elastic and viscous compliance tensors of the viscoelastic component, 1 EA is the elastic compliance tensor that is in parallel to the viscoelastic one and v is the velocity field. One should notice that a stress- velocity formulation is being used that presents several advantages over more classical stress-displacement approaches, including the ease with which dynamic effects may be considered in the analysis. Therefore, compatibility relations are written in terms of strain velocities as G G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 353   1 ( ) 2 s Tv v v v       (2) whereas the dynamic equilibrium reads div v g      (3) Truly-mixed variational formulation By observing that the total stress  may be additively decomposed as 0 1    , the continuous variational formulation of the problem at hand may be obtained by eliminating the strain tensor  in Eq. (1) and (2) testing the resulting equation by two virtual stress fields 0 1,  and the equilibrium Eq. (3) by a virtual velocity field w so as to write: find 20 1( , , ) (div, ) (div, ) ( )v H H L        such that 0 0 0 0 0 0 0 1 1 1 1 0 1 , , , div  0 , , div  0 , div  , div  , , E V E A A v A v v w w w g w                                                (4) 2 0 1(div, ), (div, ), ( ).H H w L          In more compact form Eq. (4) may be rewritten as . 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 T E V T E A A B A B M v B B v g                                                  (5) It should be emphasized that static and dynamic models differ only in that the mass matrix M may be neglected or considered in statics or dynamics, respectively. Eq. (5) amounts to a system of ordinary algebraic-differential equations (ADE) in the former case, of differential equations in the latter. Imposing Dirichlet and Neumann boundary conditions Within the adopted truly mixed formulations, Dirichlet boundary conditions on the velocity field , on uv v  (6) are imposed weakly thanks to the line integral   u n v ds    (7) that appears on the right hand side of Eq. (4)1 and (4)2, and in fact, taking no action on a boundary amounts to imposing a homogeneous condition 0, on uv   . Conversely, Neumann traction boundary conditions n t   are to be imposed strongly. However, when an additive stress decomposition is adopted, as is the case herein, i.e. 0 1    , neither 0 nor 1 are known but their sum. Therefore a Lagrange multiplier approach is adopted to enforce weakly a condition of type 0 1( ) n t    . Therefore, the resulting linear algebraic-differential system takes the form . 0 0 0 2 0 1 1 12 2 2 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 T T E V T T E A A B B A B B v v gM B B tB B                                                                 (8) G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 354 that may be compactly rewritten as 1 0D Y D Y F  (9) with an obvious definition of vectors and matrices. The Arnold-Winther finite element The triangular Arnold-Winther finite element used in this paper is the lowest-order of the family of finite elements introduced in the pioneering paper [5]. Fig. 2 shows the relevant degrees of freedom that may be listed as follows. Arnold-Winter Stress DOFS Dispacement DOFS Figure 2: Degrees-of-freedom of the Arnold-Winther stress element. As to the stresses, one should notice that the symmetry of the stress tensor is imposed strongly so that the components to be approximated are 11 22 12, ,   and one ends up with 24 degrees-of-freedom:  the three components of the stress tensor 11 22 12, ,   at each vertex of the triangle (9 dofs),  the moments of order zero and one of the traction vector n  along each edge of the triangle (12 dofs),  the averages of the components of the stress tensor over the triangle, i.e. 11 22 12, , T T T      (3 dofs). As to the displacements, the two components xu and yu are linear on each element and globally discontinuous. Implementation details No doubt that implementing the Arnold-Winther finite element represents a severe challenge, especially if the goal of minimizing the memory storage is pursued as it should if one recalls that this element is far more expensive than more conventional ( , )u p elements and other ( )H div elements such as the Johnson-Mercier element [10]. A possible implementation of the Arnold-Winther element is proposed in [6] that exploits indirect evaluations of the relevant stiffness matrices. Within the present paper, a different semi-analytical approach has been followed inspired by classical isoparametric elements according to which stress shape functions are first computed analytically on a parent triangular domain. Though not being isoparametric, the Arnold-Winther element enjoys the well-known Piola transformation property ˆ ˆ( ) ( ) Tx B x B  (10) where B is the Jacobian of the affine transformation between the reference and actual configuration. This allows to compute the stiffness matrix in any actual deformed configuration. G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 355 FORMULATION AND DISCRETIZATION (THIN BEAM CASE) he present section mimics the content of previous one but for the fact that it focuses on a Bernoulli viscoelastic beam. Truly mixed variational form and FEM discretization The phenomenological model of Fig. 1 is now used to express the adopted viscoelastic behavior in terms of bending moments M and dual curvatures  . Eq. (1) is then rewritten in the form 0 00 0 11 1 1 1 E V E M M EJ EJ M EJ           (11) and coupled to the equilibrium equation '' M q  (12) to arrive at the following Hellinger-Reissner truly-mixed formulation: find 2 2 20 1( , , ) ( )M M v H H L   such that: * * * 0 0 0 0 00 0 0 0 0 * * 1 1 10 0 0 0 1 0 0 0 0 1 1 0 1 0 E V E M M M M v M EJ EJ M M v M EJ vw M w M w qw                                              (13) * 2 * 2 2 0 1, ,M H M H w L      , which is apparently the thin beam counterpart to Eq. (4). As to the finite element discretization, the bending moment is approximated by means of cubic Hermite polynomials so that each node has two degrees of freedom, i.e. the bending moment itself and its first spatial derivative, i.e. the shear value. By analogy with the continuum case, the velocities are interpolated by elementwise linear, globally discontinuous polynomials. For later use, Eq. (14) is rewritten in compact form as 1 0D Y D Y F  (15) where 1 * * * 0 0 0 0 00 0 0 0 0 * 1 * 1 1 1 00 0 0 0 0 0 0 1 1 0 0 0 1 0 0 , 0 0 0 0 0 E V E M M M M v M EJ EJ D M M D v M EJ vw M w M w                                                                 (16) and 0 1 0 , 0 M Y M F v q                    T G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 356 EIGENVALUE-BASED TOPOLOGY OPTIMIZATION AND SOLUTION OF THE DIFFERENTIAL-ALGEBRAIC PROBLEM Eigenvalue based topology optimization he abstract unforced equation of motion for the two systems intriduced above reads 1 0 0D Y D Y  (18) that, by considering a vector solution of type ( ) stY t Ye , gives rise to the classical generalized eigenproblem 0 1[ ] 0.D sD Y  (19) The class of eigenvalue based optimal design problems considered herein may be written as 0 1 0 max ( ) s.t. [ ] 0 0 F s D sD Y V                 (20) where ( )F s is a scalar-valued objective function depending on the eigenvalues s of the problem, (20)2 is the eigenproblem that enters the optimization procedure as a constraint and (20)3 and (20)4 are global and local material density bounds, respectively. A classical SIMP model relating elastic and viscoelastic moduli to the third power of the material density is considered whereas a linear dependence between the mass density and the material density is adopted, i.e. the two densities coincide as a matter of fact. The solution of the topology optimization problem is sought by means of a sequential quadratic programming approach that has been implemented in Matlab in a self-consistent code that include the finite-element analysis part as well. Numerical solution of the Algebraic-Differential problem The solution of the (forced) algebraic-differential equation that governs the dynamics of the viscoelastic problem 1 0D Y D Y F  (22) is far from being trivial, especially in the 2D case wherein the algebraic part of the system has considerable dimension because of the need to impose several Cauchy boundary conditions on all the unconstrained sides of the domain. For this paper’s sake, a second-order accurate numerical approximation consisting of a trapezoidal rule followed by a 2-step backward difference scheme, as suggested in [4], has been adopted: 1 1 1 2 2 2 1 0 11 1 1 12 1 0 1 2 2 3 1 2 2 2 2 n n nn n n nn n n n t D Y Y f f D Y Y t D Y Y Y D Y f                                            (23) NUMERICAL STUDIES Analysis and topology optimization of a clamped/simply-supported beam clamped/simply-supported beam of length 8L  is investigated first. The initial design consists of a uniform cross-section and density. If x denotes the local density of the beam (that may be interpolated element wise or the node level), the following SIMP-like interpolations are adopted: T A G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 357 0 0 0 0 min max min 1 1 1 1 min max min 0 0 0 0 min max min min max min ( ) ( ) ( ( ) ( )) ( ) ( ) ( ( ) ( )) ( ) ( ) ( ( ) ( )) ( ) ( ) ( ( ) ( )) p E E E E p E E E E p V V V V EJ x EJ x x EJ x EJ x EJ x EJ x x EJ x EJ x EJ x EJ x x EJ x EJ x x x x x x                       (24) where 3p  and the minimum and maximum values adopted in the numerical simulations are respectively given in Tab. 1 and 2. 0 min( )EEJ x 1 min( )EEJ x 0 min( )VEJ x min( )x MIN VALUES 0.1 0.05 0.005 0.02 MAX VALUES 10 5 0.5 3 Table 1: SIMP interpolation: minimum and maximum values. Before introducing a few choices for ( )F s to be tested numerically, a few comments on the spectral properties of the problem at hand are in order. Notice in fact that neither is the system classically damped, nor enjoys the classical structure of state-space structural dynamics whose eigenproperties are well known. Reference is made to the doubly-clamped beam but similar considerations apply to all constraint cases as well. If NUMEL denotes the number of elements, the total dimension N of the problem may be written and decomposed as Velocity dofsMoment dofs Viscoleastic law dofs 2 ( 1) 2 2 ( 1)N NUMEL NUMEL NUMEL           (21) The 2 NUMEL constitutive laws (two devices in parallel per element) bring into the formulation an equal number of purely real eigenvalues, whereas moments and velocities dofs, globally amounting to 4 2NUMEL  , are associated to 2 NUMEL pairs of complex and conjugate eigenvalues and 2 null eigenvalues (that become 3 in the case of a clamped- supported beam). As to the objective function ( )F s in Eq. (20), the following three cases are considered:        1 1 1 2 2 2 2 1 2 2 3 3 3 2 2 1 ( ) Im( ) ( ) Im( ) Im(s ) ( ) Im( ) Im(s ) Im( ) Im(s ) F s s F s s F s s s              (25) where 4 6 31 2 310 , 10 , 10     . The first objective is very classical and aims to the maximization (of the imaginary part of) the lowest eigenvalue. The second objective aims to maximize the spectral band between the first two adjacent eigenvalues to avoid dangerous mode superpositions. The third function follows the same streamline as the second but involves the third eigenvalue as well. For the three optimal design problems with objectives 1 2 3, ,F F F , Fig. 3, 4 and 5 respectively present the convergence path (on the left) and the optimal density distribution along the beam axis (on the right). Fig. 6 and 7 show the first three velocity and moment eigenmodes, respectively, and demonstrate the strict correlation between the optimal density at convergence and the spectral properties of the beam. Tab. 3 shows a few values for the three design cases. The three objective functions 1 2 3, ,F F F are respectively dominated by the first, second and third eigenmode and therefore each design maximizes its own objective (boldface values on the diagonal in Tab. 3) but, with no surprise, happen to lead to small values when tested on a different objective than the one used in the optimization (first three columns in Tab. 3). G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 358 0 200 400 600 800 1000 1200 1400 1600 1800 −200 −180 −160 −140 −120 −100 −80 −60 −40 −20 Number of Iterations O b je ct iv e F u n ct io n F 1 (s ) 0 1 2 3 4 5 6 7 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Beam axis O p tim a l M a ss D e n si ty Figure 3: Objective function 1: path to convergence and optimal density. 0 500 1000 1500 2000 2500 −2500 −2000 −1500 −1000 −500 0 Number of Iterations O b je ct iv e F u n ct io n F 2 (s ) 0 1 2 3 4 5 6 7 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Beam axis O p tim a l M a ss D e n si ty Figure 4: Objective function 2: path to convergence and optimal density. 0 500 1000 1500 2000 2500 −7000 −6000 −5000 −4000 −3000 −2000 −1000 0 Number of Iterations O b je ct iv e F u n ct io n F 3 (s ) 0 1 2 3 4 5 6 7 8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Beam axis O p tim a l M a ss D e n si ty Figure 5: Objective function 3: path to convergence and optimal density. G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 359 0 1 2 3 4 5 6 7 8 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Beam axis D is p la ce m e n t E ig e n m o d e 1 0 1 2 3 4 5 6 7 8 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Beam axis D is p la ce m e n t E ig e n m o d e 2 0 1 2 3 4 5 6 7 8 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 Beam axis D is p la ce m e n t E ig e n m o d e 3 Figure 6: First three velocity eigenmodes. 0 1 2 3 4 5 6 7 8 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 x 10 −5 Beam axis B e n d in g M o m e n t E ig e n m o d e I 0 1 2 3 4 5 6 7 8 −12 −10 −8 −6 −4 −2 0 2 x 10 −5 Beam axis B e n d in g M o m e n t E ig e n m o d e I I 0 1 2 3 4 5 6 7 8 −6 −4 −2 0 2 4 6 8 10 x 10 −5 Beam axis B e n d in g M o m e n t E ig e n m o d e I II Figure 7: First three moment eigenmodes. DESIGN OBJECTIVE 1 ( )F s 2 ( )F s 3( )F s 1Im( )s 2Im( )s 3Im( )s 1( )F s 0.019 0.198 1.67 0.019 0.034 0.072 2 ( )F s 0.0065 2.432 2.945 0.0065 0.0558 0.0785 3( )F s 0.0055 0.189 6.261 0.0055 0.0192 0.097 Table 3: Numerical results. To assess accuracy and convergence of the proposed approach reference is made to Fig. 8, 9 and 10 that, for the objectives 1 2 3, ,F F F respectively, show the optimal values taken on by the objective functions versus the number of elements used for the discretization (on the left) and relevant optimal designs (on the right). The method is clearly robust in that optimal solutions are clearly mesh insensitive and even using a coarse mesh of 8 elements a good approximation to the exact optimal solution (that is of course analytically unknown and is herein assumed to coincide with the numerical one computed with a mesh of 64 elements) is found. Of course, the higher the mode number that dominates the objective function the finer the mesh needed to find an accurate solution: one may in fact see that the solution using 8 elements (Fig. 8 on the right) nearly coincides with the exact one when the objective function is 1F that depends on the first eigenmode only, whereas the accuracy decreases when solving the optimal design problems governed by 2F and 3F that involve higher eigenmodes. However, this should not be seen as a limitation of the proposed optimization approach in that when a coarse mesh is not capable to approximate accurately the eigenmodes for a given material density, the solution of the optimization procedure cannot be any better. G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 360 0 10 20 30 40 50 60 70 0.0165 0.017 0.0175 0.018 0.0185 0.019 0.0195 0.02 Number of Elements F 1 (s ) 0 1 2 3 4 5 6 7 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Beam axis O p tim a l D e n si ty − F 1 (s ) 23 elements 24 elements 25 elements 26 elements Figure 8: 1 sF : Convergence vs number of elements (left) – Optimal densities vs number of elements (right). 0 10 20 30 40 50 60 70 −2.45 −2.4 −2.35 −2.3 −2.25 −2.2 −2.15 −2.1 −2.05 −2 −1.95 Number of Elements F 2 (s ) 0 1 2 3 4 5 6 7 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Beam axis O p tim a l D e n si ty − F 2 (s ) 23 elements 24 elements 25 elements 26 elements Figure 9: 2 sF : Convergence vs number of elements (left) – Optimal densities vs number of elements (right). 0 10 20 30 40 50 60 70 −7 −6 −5 −4 −3 −2 −1 Number of Elements F 3 (s ) 0 1 2 3 4 5 6 7 8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Beam axis O p tim a l D e n si ty − F 3 (s ) 23 elements 24 elements 25 elements 26 elements Figure 10: 3 sF : Convergence vs number of elements (left) – Optimal densities vs number of elements (right). G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 361 Analysis of a plane cantilever The plane cantilever of Fig. 11, uniformly loaded on the right side is considered. The spatially uniform load is modulated in time by means of the function 1 1 sin 2 ( ) 0 t t t tF t else     (26) where 1 2t s and [0,10] .t s As to the physical properties of the structure the following values are adopted as to the isotropic elastic and viscous compliance tensors that are defined in terms of Lamé constants: 0 1 0 100 80 100 75 90 30 E E V A A A         (27) A uniform mass density 0.5  is further considered. Fig. 12, 13 and 14 show the maps of the stress components , ,xx xy yy   at final time 10t s . Figure 11: Plane viscoelastic cantilever under investigation. Figure 12: xx (normal stress). G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 362 Figure 13: xy (shear stress) Figure 14: yy (vertical stress) To appreciate the diffusion of the stress components with time, Fig. 15 and 16 show a few snapshots that display the variation of the stress component ,xx xy  at regularly spaced time stations. Figure 15: Evolution of xx (normal stress) with time. Figure 16: Evolution of xy (shear stress) with time. Finally, Fig. 17 shows the time-variation of the stress xx at the upper-left corner of the cantilever. CONCLUSIONS truly-mixed variational formulation for the analysis of viscoelastic continua and thin beams has been presented that uses stresses (moments) as primary variables and velocities (instead of displacements) as Lagrange multipliers. A topology optimization method for viscoelastic devices has then been developed and applied to thin beams within a few representative numerical simulations. Ongoing extensions include applications to two dimensional viscoelastic devices with respect to eigenvalue-based design objectives. A G. Maurelli et alii, Frattura ed Integrità Strutturale, 29 (2014) 351-363; DOI: 10.3221/IGF-ESIS.29.31 363 Figure 17: Evolution of xx (normal stress) with time at the upper left corner of the cantilever. REFERENCES [1] Phan-Thien, N., Understanding viscoelasticity, Advanced Texts in Physics, Springer-Verlag, Berlin, (2002). [2] Bruggi, M., Venini, P., Eigenvalue-based optimization of incompressible media using mixed finite elements with application to isolation devices, Computer Methods in Applied Mechanics and Engineering, 197(13-16) (2008) 1262- 1279. [3] Brezzi, F., Fortin, M., Mixed and hybrid finite element methods, Springer Series in Computational Mathematics. Springer-Verlag, New York, 15 (1991). [4] Rognes, M.E., Winther, R., Mixed finite element methods for linear viscoelasticity using weak symmetry', Math. Models & Methods in Applied Sciences, 20(6) (2010) 955-985. [5] Arnold, D., Winther, R., Mixed finite elements for elasticity, Numer. Math., 92 (2002) 401-419. [6] Carstensen, C., Gunther, D., Reininghaus, J., Thiele, J., The Arnold-Winther mixed FEM in linear elasticity. Part I: Implementation and numerical verification, Computer Methods in Applied Mechanics and Engineering, 197 (2008) 3014-3023. [7] Diaz, A., Kikuchi, N., Solution to shape and topology eigenvalue optimization problems using a homogenization method, Int. J. Numer. Methods Engrg. 35(3) (1992) 1487–1502. [8] Krog, L.A., Olhoff, N., Topology optimization of plate and shell structures with multiple eigenfrequencies, In: N. Olhoff, G.I.N. Rozvany, Proceed. WCSMO-1 First World Congress of Structural and Multidisciplinary Optimization, 675–682, (1995). [9] Pedersen, N.L., Maximization of eigenvalue using topology optimization, Struct. Multidiscip. Optim. 20 (2000) 2–11. [10] Johnson C., Mercier, B., Some equilibrium finite elements methods for two dimensional elasticity problems, Numer. Math., 30 (1978) 103–116. [11] Bruggi, M., Venini, P., A mixed FEM approach to stress-constrained topology optimization, Int. J. Numer. Methods Engrg., 73 (2008) 1693–1714.