An intuitive computational multi-scale methodology and tool for the dynamic modelling of viscoelastic composites and structures

Composite Structures 144 (2016) 131–137

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

An intuitive computational multi-scale methodology and tool for the dynamic modelling of viscoelastic composites and structures M. El Hachemi a,⇑, Y. Koutsawa a, H. Nasser a, G. Giunta a, A. Daouadji c, E.M. Daya b, S. Belouettar a a

Luxembourg Institute of Science and Technology, LIST, 5, Avenue des Hauts Fourneaux, L-4362 Esch/Alzette, Luxembourg Laboratoire d’Etude des Microstructures et de Mécanique des Matériaux, LEM3, UMR CNRS 7239, Université de Lorraine, Ile du Saulcy, F-57045 Metz, France c Université de Lyon-INSA, Laboratoire de Génie Civil (LGCIE), Batiment Coulomb, F-69100 Villeurbanne, France b

a r t i c l e

i n f o

Article history: Available online 20 February 2016 Keywords: Multi-scale method Viscoelastic Homogenisation Vibration Composites

a b s t r a c t This paper proposes an intuitive computational multi-scale homogenisation procedure and tool for the estimation of the effective static and mechanical properties of complex viscoelastic composite material and structures. The proposed solution consists in computing numerically the complex effective properties (storage and loss moduli) as function of frequency. The developed numerical tool is coupled with ABAQUS FEA in order to ease the uptake of the technology code and allows for an accurate and automatic simulation of composite engineering structures with substantially less human intervention and a rational control of the error. Details concerning the numerical implementation and a selection of representative numerical examples are provided. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Structural vibrations control is of primary importance for enhancing safety and improving structures and systems performances. The most prevalent strategy for structural vibration damping is the one that consider adding viscoelastic inhomogeneities either as thin layers or as dilute inclusions depending on the practical considerations. These composite materials and structures yield superior vibrational energy absorption capability. They particularly offer the advantage of high damping, low cost, ease of implementation with low weight. Analytical studies were devoted to simple structures and finite element simulations were introduced to design structures with complex geometries and generic boundary conditions. Most analyses use a complex frequency-independent dynamic modulus to describe the rheological behaviour of viscoelastic materials. However, it is widely known from engineering practice, that the storage and loss moduli of viscoelastic materials are strongly frequency and temperature dependent. This dependence is especially significant for elastomers and polymers, such as those used in laminated composites. A comprehensive review of the various research methods and theory calculation models could be found in the very recent review by Zhou [1].

⇑ Corresponding author. Tel.: +352 275 888 1; fax: +352 275 885. E-mail address: [email protected] (M. El Hachemi). http://dx.doi.org/10.1016/j.compstruct.2016.02.032 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

Nevertheless, viscoelasticity is not the only mechanism for damping. As a matter of fact, defects and interfaces can also contribute to damping vibrational dissipating energy. Therefore, the micro-structure greatly affects the damping ability of a material. Consequently, new materials with tailored micro-structure have attracted increased interest because of the possibility to control both their stiffness and damping characteristics i.e. loss modulus. In parallel, the recent advancements in micron-level additive manufacturing have enabled unprecedented accuracy in controlling the material architecture resulting in an unrestricted number of possibilities for multi-scale material design. In view of acceptance of these new composite materials for vibration damping, the development of models and numerical solutions and tools is obvious. It is through these models that the effect of shapes, sizes and location of phase inclusions on the dynamical properties can be investigated. These models should be capable of correlating the micro-structural response with the overall macroscopic dynamical behaviour since the applied loads are at the structural level. A series of micro-mechanical approaches for the estimation of the damping properties of composite materials, as a function of the properties of the constituent materials have been developed in [2–12]. Generally, the effective behaviour of the viscoelastic composite in the time domain is obtained using a combination of the correspondence principle with the homogenization solution of the corresponding elastic problem. Then, the solution of the viscoelastic effective properties is obtained using the Laplace Carson transform technique. This permits transferring a viscoelastic

132

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

problem into a symbolic elastic problem with known solution. The main technical issue when using these analytical methods towards up-scaling of viscoelastic properties of composites, is performing the inverse Laplace transform since the expressions of the effective properties in the Laplace space are not rational fractions as functions of the Laplace variable. Even sophisticated approaches have been developed and used, such as in [13,14], to perform the inverse transform, basically, all these analytical and semi analytical methods are restricted to canonical micro-structure shapes like sphere or cylinder like shapes. Recent works proposed the use finite element method (FEM) solution in multi-scale context where the macroscopic parameters are obtained by volume averaging over statistically representative volume which are then used at the macroscopic level [15–21]. This latter numerical framework can be an alternative approach to answer the previous drawbacks. Typically, this approach is based on a hierarchical decomposition of the solution space into a local solution and a global one and by enforcement of the solution compatibility conditions. A multilevel finite element methodology is introduced where the hierarchical character of model description and simulation results are exploited to expedite the analysis of problems. This makes it ideal for local/global analyses where solutions from a local model, i.e. the micro-structure, are used to derive the solution for the global model of the structure. In this paper, an intuitive numerical procedure for the estimation of the effective properties of the viscoelastic composites is developed and a computational multilevel methodology and tool for the prediction of dynamical properties of composite structures is proposed as an alternative to the direct simulation which requires enormous computing resources. The proposed solution consists of computing numerically the complex effective proprieties (storage and loss moduli) as function of frequency. This is obtained using a steady state dynamic direct solver with a sweep over the frequency band of interest and then simulating the Mechanical Impedance (MI) [22] test for an equivalent homogeneous structure in order to predict the spectral response in terms of driving point mobility. When compared to the direct numerical simulation, this method ensures a considerable saving in computational effort. The method is suitable to deal with composites with complex micro-structure and to cope with the inaccuracy of traditional constitutive relations. It is also suitable and effective to simulate composites with high volume fraction. Hence, it allows for tailoring the effective properties to specific application requirement. The developed numerical tool is coupled with ABAQUS FEA in order to ease the uptake of the technology code and allow for automatic simulation of complex 3D composite engineering structures and boundary conditions. Details concerning the numerical implementation and a selection of representative numerical examples are provided. 2. Multi-level FEM approach and problem formulation The multilevel principle assumes that if a heterogeneous structure ðXÞ contains evenly distributed micro-topologies, a Representative Volume Element (RVE) ðXRVE  XÞ can be, then, defined where the averaged micro-structural behaviour can be smeared as an homogeneous one over the macroscopic domain ðXÞ (see Fig. 1). A concise description of the particular technique of multilevel finite element method, that constitutes the first significant part of the approach described in the current work, can be found in [17,19,21]. 2.1. Problem at the macro-scale The problem at the macro-scale is governed by the equation of motion, which is an expression of Newton’s second law:

Fig. 1. Schematic representation of the multilevel multi-scale approach.

q


tÞ @ 2 Uðx;
tÞ
ðx; tÞ ¼ Fðx; þrr @t2

ð1Þ


¼ ½u
1 ; u
2 ; u
3  is the displacement vector, the bar sign ð
Þ where U denotes the macroscopic value, x is the vector position and t is the time. q stands for density. The composite material is considered as viscoelastic, thus the stress rðx; tÞ is related to the strain ðx; tÞ by the Boltzmann convolution integral:

r
ðx; tÞ ¼

Z 0

t


 sÞ : @
eðx; sÞds Cðt @s

ð2Þ

while t is the instant of loading, s spans the history time and C is the relaxation modulus. If we apply the Laplace–Carson transformation, Rt  LC ½f ðtÞ ¼ f ðpÞ ¼ p 0 f ðtÞept dt, to the convolution integral in Eq. (2), the constitutive relation becomes a simple tensor contraction in the Laplace–Carson space:


ðx; tÞ ¼ r
 ðx; pÞ ¼ C
 ðpÞ : e
 ðx; pÞ LC ½r

ð3Þ


xÞ ¼ C
0 ðxÞ þ iC
00 ðxÞ Cð

ð4Þ

pffiffiffiffiffiffiffi with p ¼ a þ ix and i ¼ 1. a is an arbitrary real value and the Laplace–Carson domain is, then, the frequency domain.
xÞ is the complex fourth order tangent modulus tensor:
 ðpÞ  Cð C

where x is the angular frequency. By considering the excitation force to be a monochromatic plane wave, the displacement field in the medium is also a monochromatic plane wave:


tÞ ¼ UðxÞ
Uðx; expðixtÞ

ð5Þ

By replacing the displacement form in Eq. (5) and the constitutive law in Eq. (4) in the governing equation, Eq. (1), the time-free wave equation of elasto-dynamics is obtained:

qx2 u
i þ @ j C
ijkl @ l u
k ¼
f i

ð6Þ

This system can be solved using a numerical perturbation procedure, where the perturbed solution is obtained by a linearisation in the current base state. Structural and viscous damping can be included in the numerical procedure using the Rayleigh and structural damping coefficients. Besides, global damping coefficients can be specified at the procedure level to define additional viscous and structural damping contributions. In the context of a multilevel finite element analysis, the stress field at the macro-scale level can be computed by solving a local nonlinear finite element problem.

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

133

2.2. Problem at the micro-scale

3. Resolution strategy and computation procedure

The material is considered heterogeneous at the micro-scale with an evenly distributed inclusions so that the micro-structure can be characterised by a unit cell that occupies a domain XRVE . The static and dynamic equilibrium equations are given as:

The implemented two scales homogenisation strategy is based on the formulation described above, see Fig. 1. The flowchart of the numerical procedure including the coupling to ABAQUS FEA is reported in Fig. 2. A micro-structure simulation tool is developed to automatically generate various 2D and 3D micro-structure models within a prescribed geometrical domain. The size of the geometrical domain should be set sufficiently small compared to the operating wavelength. In addition, it should comply with the volume fraction of each of the material components. The specified domain defines, therefore, an RVE. The used mesh size is fine enough to represent accurately all the geometrical features of the micro-structure. To obtain the homogenized effective composite material properties, the periodic boundary conditions in Eq. (8) are applied to the RVE by coupling opposite nodes on the opposite boundary surfaces. In order to apply these boundary conditions on the Finite Element (FE) model of the RVE, the meshes on the opposite surface must be same. To do so, a micro-domain tessellation has been generated such that each pair of triangles meshes of opposite faces of the RVE are identical images. A constraint equation is imposed for each displacement component at two corresponding nodes. The tool generates all the required constraint equations. This provides an efficient and powerful approach to construct a FE model for composite materials simulation with a great variety of micro-structures. In order to derive the long
1 , a static linear perturbation solver is used. term elasticity tensor C The needed averaged static stress tensor is evaluated as outlined in
xÞ, a steady (11). In order to obtain the dynamic elasticity tensor Cð state direct linear perturbation solver is used. Similarly, the needed averaged dynamic stress tensor is evaluated using the same Eq. (11) with sweep over the frequency domain of interest. It should be stated that for an isotropic material, the effective properties, hence obtained, are used for solving the system of Eq. (6) in the frequency domain using steady state direct linear perturbation solver with a sweep over the frequency range of interest in order to calculate the acceleration at the drive point AðxÞ.



rrðxÞ ¼ 0 : for static loading; qx2 U þ rrðxÞ ¼ 0 : for harmonic loading

ð7Þ

where rðxÞ is the Cauchy stress tensor at a micro-scale point x. The total strain is computed in ABAQUS FEA code by approximating the integral of the strain rate over an increment using a central difference algorithm. Periodic boundary conditions are adopted on the limits of the RVE such as:


ðxþ  x Þ
 ¼ E
þ  u u

ð8Þ


 represent the imposed displacement at the oppoIn Eq. (8), u site boundaries of the RVE. The exponents ðÞ are associated with node indexes on opposite sides of the RVE, x are the nodes posi
is the imposed deformation. Small deformations are tions and E considered and a linear constitutive relation is assumed for each phase:

rðx; xÞ ¼ C ðiÞ ðxÞ : eðxÞ;

ð9Þ

where C ðiÞ refers to the fourth-order elastic tensor associated with phase (i) and e is the Green–Lagrange strain tensor. The weak form associated with the micro-scale problem in Eq. (7) is written as: Find U that satisfies Eq. (8) such that:

(R

R

R

r : rx du dX ¼ XRVE r : de dX ¼ 0 R R 2 XRVE r : rx du dX ¼ XRVE r : de dX ¼  XRVE qx u  dudX XRVE

ð10Þ

d stands for a virtual variation. The effective stress is evaluated as field average of the local response through the following relation:

hrðxÞi ¼

1 V

Z

XRVE

rðxÞdX

ð11Þ

where V denotes the volume of the RVE. As reported in Özdemir et al. [23], the usual relationships between stress measures (e.g.
are not the Cauchy and the first Piola–Kirchhoff stress tensors P) always valid for the volume averages of the micro-structural coun
the Cauchy stress tensor on terparts. If the averaging is based on P, the macro-scale should be defined as:

hrðxÞi ¼

1

 P  F det F


4. Integrated micro-structure generation Due to the recurrent use of the homogenisation process, an automatic procedure has been coupled to ABAQUS FEA code as shown in Fig. 3. A similar approach has been adopted in Yuan and Fish [19] and a comprehensive procedure on the implementation of such

ð12Þ

where F
is the deformation gradient tensor. In ABAQUS, the averag
at the micro-scale is performed after ing of the Kirchhoff tensor P converting the obtained micro-scale Cauchy tensor into the first Piola–Kirchhoff tensor. It should be noticed that when the geometric non-linearity is not considered, the deformation gradient tensor F
at the micro-scale is equal to the identity tensor I and, thus, the Cauchy stress tensor and the first Piola–Kirchhoff tensor are equal. Hence, it is not necessary to use Eq. (12) to obtain the ‘true’ value of the Cauchy stress tensor. The averaging defined in Eq. (11) is sufficient in this case to obtain the homogenized constitutive relations. Eq. (11) can be rewritten as:


hrðx; xÞi ¼ C
eff : E

ð13Þ

The two problems at the micro and macro scales are coupled through the following relation that enables the evaluation of the effective elastic tensor:

@ rðx; xÞiij
eff ¼ h C ijkl
kl @E

ð14Þ Fig. 2. Flowchart of the proposed multi-scale simulation framework.

134

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

the complexity of the RVE micro-structure, these two conditions can be difficult to full-fill, in this case, a locally refined mesh should be used. At the first calculation step, the boundary condition for each loading step fT 11 ; T 22 ; T 33 ; S12 ; S13 ; S23 g are applied using three dummy nodes with three degrees of freedom each:

h i udj i ¼ dij ETii þ ð1  dij ÞESij DX i

ð15Þ

where ETii and ESij correspond to the prescribed traction (T) and shear th

(S) loading steps. DX i is the length of the i dimension. udj i denotes the ith displacement component of the jth dummy node. Afterwards, a linear multi-point constraint is used to impose the periodic boundary conditions. Three different sets of nodes are identified: nodes on opposite faces (‘f ’) in the traction direction T jj : fþ

f

ui j  ui j þ uidj ¼ 0 i; j ¼ 1; 2; 3

ð16Þ

nodes on diametrically opposite edges (‘e’) in the shear direction Sij : eþ

e

ui ij  ui ij þ uidi þ uidj ¼ 0 i; j ¼ 1; 2; 3

ð17Þ

nodes on diametrically opposite corners (‘c’) in the shear direction Sij ; Sik Fig. 3. Sketch of the RVE micro-structure generation tool.

computational multilevel approach in ABAQUS is proposed in Tchalla et al. [21]. The first step of the procedure involves the creation of a valid geometric representation of the micro-structure, see Figs. 3 and 4. The considered shapes are mainly basic geometries: brick, tetrahedron, cylinder, ellipsoid and sphere. The heterogeneities can be modelled as a random distribution with respect to their content and shape inside the material as depicted in Fig. 4. The resulting geometry and materials are used as inputs to the mesh generator module of ABAQUS. As previously mentioned, a suitable meshing tool that produces imaged triangulation of the opposite faces of the RVE is developed in order to apply the required periodic boundary conditions in Eq. (8). The quality of tetrahedra elements should satisfy some constraints. For instance, the SteadyStateDynamicDircet solver requires that the shape factor (the ratio between a tetrahedron volume and the volume of the equilateral tetrahedron with the same circumsphere) and the aspect ratio (the ratio between the smallest and longest edges of the same tetrahedron) should be higher than 0.1. Depending on

cþ

c

ui iðj;kÞ  ui iðj;kÞ þ uidi þ uidj þ uidk ¼ 0 i; j; k ¼ 1; 2; 3

ð18Þ

The instruction lines are written in an input file and loaded into ABAQUS. In the post-processing step, the tangent matrix in Eq. (14) is computed by solving a system of six equations corresponding to the six loading steps. The effective proprieties are then obtained in ABAQUS format as shown in Fig. 5b. 5. Numerical application 5.1. Case of spherical viscoelastic inclusion The numerical homogenisation tool is used to estimate the effective loss factor of the composite material. The composite matrix is made of purely elastic material (E ¼ 2:14 GPa; m ¼ 0:3 and g ¼ 0). The spherical inclusion made of a viscoelastic Epoxy resin. The dynamic Young modulus and loss factor versus frequency of the Epoxy are shown in Fig. 6. The obtained effective loss factor is compared with results obtained using a linear averaging approach:

gav g ðxÞ ¼ f ginc ðxÞ þ ð1  f Þgmat

ð19Þ

where f is the volume fraction of the inclusion. For this calculation, three volume fractions were considered: 10%; 20% and 30%. Results are shown in Fig. 7. The difference is about two orders of magnitude between the values (FEA) computed using the developed computational methodology and the analytical one (Avg). Therefore, it is important to state that the estimation of the effective loss factor of the composite using averaging formulas can yield inaccurate results. 5.2. Effect of the stiffness ratio on the effective dynamical parameters In order to increase the effective loss factor of the composite, the optimal stiffness ratio between matrix and inclusion is investigated. RVE simulation tests are performed by varying the long term

Fig. 4. Example of an RVE obtained using the developed micro-structure generation tool.

Young modulus of the inclusion Einc ¼ ½106 ; 10þ6 Emat while the loss factor is kept the same as the viscoelastic Epoxy material. The volume fraction is equal to 20%. Fig. 8 shows the cross-cut

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

135

Fig. 5. (a) Cross cut contours of the max-principal stress from the RVE analysis and (b) resulting effective parameters.

Fig. 6. Viscoelastic properties of the Epoxy resin inclusions.

Fig. 8. Cross-cut contours of shear stress S11 from the FE analysis ðEinc =Emat  1Þ.

effective loss factor can be maximised if the stiffness of the matrix and inclusion are equal. 5.3. Hollow spherical viscoelastic inclusion

Fig. 7. Effective loss factor for spherical viscoelastic inclusion: comparison between analytical (Avg) and numerical homogenisation (FEA).

contours of the shear stress S11 from the FE analysis. Fig. 9 presents the variation of the effective loss factor versus the stiffness ratio between the matrix and inclusion. Results clearly indicate the

The loss factor, which describes the energy dissipation ability, should be high for an effective vibration damping. The precedent results show that the effective loss factor depends on the dynamic response of the micro-structure rather than the viscoelastic content. This fact is further investigated by considering an hollow Epoxy sphere as shown in Fig. 10a. Details of the cross-cut contours of the max-principal stress are shown in Fig. 10b. The inclusion has a core of radius ra made of the same material as the matrix and a coating made of viscoelastic Epoxy layer of thickness r b  r a . Four thickness values are considered, ra ¼ 0; 14 ; 12 and 34. Fig. 11 presents the effect of the viscoelastic r b

thickness on the dynamic modulus. The trend is in-line with the

136

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

10,0000%

1,0000%

Loss Factor [%]

0,1000%

0,0100%

f=100 Hz f=300 Hz f=500 Hz f=700 Hz f=900 Hz f=1100 Hz f=1500 Hz f=2000 Hz

0,0010%

0,0001%

0,0000% 1,E-06

1,E-04

1,E-02

1,E+00 1,E+02 Rao Einc/Emat

1,E+04

1,E+06 Fig. 11. Comparison between different thickness of the viscoelastic hollow sphere and effect on dynamic moduli.

Fig. 9. Evolution of the effective loss factor at different frequencies when varying stiffness ratio between the matrix and inclusion.

Epoxy content. As the Epoxy is relatively soft compared to the elastic matrix, the lower the Epoxy content, the higher the dynamic modulus. The curves that show the effect of the Epoxy layer thickness on the loss factor are presented in Fig. 12. The trend is inverse to the viscoelastic Epoxy content: rra ¼ 34 exhibits the highest effecb

tive loss factor. 5.4. Two-scale vibration analysis of heterogeneous strip Finally, the two-scale simulation approach is validated by comparing the results from the simulation of the vibration of a heterogeneous rectangular strip with experimental measurements. The strip has the dimension of 300 30 13:2 mm and consists of composite mixture, with pure elastic matrix of properties: Young’s modulus E ¼ 2:14 GPa, mass density q ¼ 696:9 kg=m3 and Poisson’s ratio m ¼ 0:3. The inclusions are viscoelastic with the proprieties shown in Fig. 6, represented by spheres of diameter 2 mm and uniformly distributed within the matrix with volumetric fraction of 20%. An Harmonic concentrated load is applied at the bottom centre of the strip to represent the electrodynamic shaker. The frequency range is ½100 Hz;2 KHz, which includes the two first transversal modes at 262 Hz and 1379 Hz. The pick-up point is at the driving point. Fig. 13 shows the driving point mobility which is the velocity in Fourier space normalised to the applied force amplitude.

(a)

Fig. 12. Effect of the Epoxy layer thickness on the effective loss factor of the viscoelastic hollow sphere.

The results computed from the two-scale FEM show that the frequency response of the structure compares very well with that obtained from the experimental measurements. The predicted eigen frequencies much very well the experiment albeit the

(b)

Fig. 10. Hollow spherical viscoelastic inclusion, (a) sketch showing details of the geometry, (b) cross-cut contours of max-principal stress.

M. El Hachemi et al. / Composite Structures 144 (2016) 131–137

Fig. 13. Comparison between experiment and simulation.

anti-resonance frequencies are slightly shifted due its strong dependence to pick-up point position. This example shows that the implemented method is accurate, flexible (the micromechanical model is managed by the tool), it can be easily linked to ABAQUS and it allows a considerable save in computational time when compared with full direct simulation. 6. Conclusion An intuitive homogenisation procedure and tool for frequency dependant viscoelastic composites is proposed. The estimation of the effective properties of complex viscoelastic composite material is performed using the developed computational tool. The proposed computational multi-scale methodology has shown to be effective and provides accurate information about the effect of the micro-structure of the global dynamic behaviour of composites. Numerical results have been validated towards analytical solutions as well as experimental results showing that the proposed tool is accurate and it can deal with composites with different micro-structures morphologies. The developed numerical tool is coupled with ABAQUS FEA via plug-in interface, to allow for automatic simulation of complex 3D composite engineering structures and boundary conditions. Finally, the developed tool allows for the design and tailoring of the effective properties for specific application requirements. Acknowledgement Dr. H. Nasser acknowledges the financial support of the Fonds National de la Recherche (FNR) of Luxembourg under the framework of OPTIPIEZO CORE project C11/MS/122837. References [1] Zhou XQ, Yu DY, Shao XY, Zhang SQ, Wang S. Research and applications of viscoelastic vibration damping materials: a review. Compos Struct 2016;136:460–80. . [2] Christensen R. Viscoelastic properties of heterogeneous media. J Mech Phys Solids 1969;17:23–41. http://dx.doi.org/10.1016/0022-5096(69)90011-8. [3] Wang YM, Weng GJ. The influence of inclusion shape on the overall viscoelastic behavior of composites. J Appl Mech 1992;59(3):510–8. http://dx.doi.org/ 10.1115/1.2893753.

137

[4] Koishi M, Shiratori M, Miyoshi T, Kabe K. Homogenization method for dynamic viscoelastic analysis of composite materials. JSME Int J Ser A, Solid Mech Mater Eng 1997;40(3):306–12 [anglais. . [5] Brinson L, Lin W. Comparison of micromechanics methods for effective properties of multiphase viscoelastic composites. Compos Struct 1998;41(3– 4):353–67. http://dx.doi.org/10.1016/S0263-8223(98)00019-1. . [6] Yi Y-M, Park S-H, Youn S-K. Asymptotic homogenization of viscoelastic composites with periodic microstructures. Int J Solids Struct 1998;35 (17):2039–55. http://dx.doi.org/10.1016/S0020-7683(97)00166-2. . [7] Yi Y-M, Park S-H, Youn S-K. Design of microstructures of viscoelastic composites for optimal damping characteristics. Int J Solids Struct 2000;37 (35):4791–810. http://dx.doi.org/10.1016/S0020-7683(99)00181-X. . [8] Miehe C, Göktepe S. A micro-macro approach to rubber-like materials. Part II: The micro-sphere model of finite rubber viscoelasticity. J Mech Phys Solids 2005;53(10):2231–58. http://dx.doi.org/10.1016/j.jmps.2005.04.006. . [9] Koutsawa Y, Cherkaoui M, Daya E. Multi-coating inhomogeneities problem for effective viscoelastic properties of particulate composite materials. J Eng Mater Technol 2009;131(2):021012.1–021012.11. http://dx.doi.org/10.1115/ 1.3086336. [10] Koutsawa Y, Azoti WL, Belouettar S, Martin R, Barkanov E. Loss behavior of viscoelastic sandwich structures: a statistical-continuum multi-scale approach. Compos Struct 2012;94(4):1391–7. http://dx.doi.org/10.1016/ j.compstruct.2011.11.003. [11] Azoti W, Bonfoh N, Koutsawa Y, Belouettar S, Lipinski P. Influence of auxeticity of reinforcements on the overall properties of viscoelastic composite materials. Mech Mater 2013;61:28–38. http://dx.doi.org/10.1016/j. mechmat.2013.02.002. [12] Mandel J. Cours de mécanique des milieux continus. Paris: Gauthier-Villars; 1966. [13] Selivanov MF, Chernoivan YA. A combined approach of the laplace transform and Pade approximation solving viscoelasticity problems. Int J Solids Struct 2007;44(1):66–76. http://dx.doi.org/10.1016/j.ijsolstr.2006.04.012. . [14] Rekik A, Brenner R. Optimization of the collocation inversion method for the linear viscoelastic homogenization. Mech Res Commun 2011;38(4):305–8. http://dx.doi. org/10.1016/j.mechrescom.2011.04.003. . [15] Hughes TJ. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput Meth Appl Mech Eng 1995;127(1–4):387–401. http://dx.doi. org/10.1016/0045-7825(95)00844-9. . [16] Zohdi TI, Oden J, Rodin GJ. Hierarchical modeling of heterogeneous bodies. Comput Meth Appl Mech Eng 1996;138(1–4):273–98. http://dx.doi.org/ 10.1016/S0045-7825(96)01106-1. . [17] Feyel F, Calloch S, Marquis D, Cailletaud G. FE computation of a triaxial specimen using a polycrystalline model. Comput Mater Sci 1997;9(1– 2):141–57 [Selected papers of the Sixth International Workshop on Computational Mechanics of Materials. doi 10.1016/S0927-0256(97)00070-0. ]. [18] Fish J, Shek K, Pandheeradi M, Shephard MS. Computational plasticity for composite structures based on mathematical homogenization: theory and practice. Comput Meth Appl Mech Eng 1997;148(1–2):53–73. http://dx.doi. org/10.1016/S0045-7825(97)00030-3. . [19] Yuan Z, Fish J. Toward realization of computational homogenization in practice. Int J Numer Meth Eng 2008;73(3):361–80. http://dx.doi.org/ 10.1002/nme.2074. [20] Attipou K, Nezamabadi S, Daya EM, Zahrouni H. A multiscale approach for the vibration analysis of heterogeneous materials: application to passive damping. J Sound Vib 2013;332(4):725–39. http://dx.doi.org/10.1016/j.jsv.2012.10.020. . [21] Tchalla A, Belouettar S, Makradi A, Zahrouni H. An abaqus toolbox for multiscale finite element computation. Compos Part B Eng 2013;52:323–33. http://dx.doi. org/10.1016/j.compositesb.2013.04.028. . [22] Allan LP, Piersol G, Thomas. Harris’s shock and vibration handbook. Mechanical impedance/mobility. McGraw Hill Professional; 2010. Access Engineering. doi: ISBN: 9780071508193. . [23] Özdemir I. Fe2 computational homogenization for the thermo-mechanical analysis of heterogeneous solids. Comput Meth Appl Mech Eng 2008;198(3– 4):602.