# Abstract

Most of the existing materials around us can be considered composite materials, since they are composed by several phases or components at certain spatial scale. The physical and chemical properties of composites, as occurs with structures composed by two or more materials, is defined by the response provided by their constituents. Therefore, a good characterization of the composite requires considering the performance of its components. In the last decades, several methods have been proposed with this approach to characterize composite materials, most of them based on multiscale techniques.

Nowadays, multiscale homogenization analysis is a popular topic in the simulation of composite materials. This is because the complexity of new composites demands of advanced analysis techniques for their correct characterization, and thanks to the continuous increase of computational capacity. However, the computational cost when multiscale procedures are taken to the non-linear range and are applied to real-size structures is still excessively high. In this context, this monograph presents a comprehensive homogenization formulation for an efficient non-linear multiscale modeling of composite structures.

The development of a composite multiscale constitutive model is addressed from two different homogenization approaches. The first one corresponds to a phenomenological homogenization procedure for the non-linear analysis of carbon nanotubes reinforced composites. The second one is a general two-scale homogenization procedure to analyze three-dimensional composite structures.

Carbon nanotubes (CNTs) have been regarded as ideal reinforcements for high-performance composites. The formulation developed takes into account explicitly the performance of the interface between the matrix and the CNTs. The load is transferred to the nanotubes through the considered interface. The composite non-linear behavior results from the non-linearities of its constituents, and in case of interface damage, it also becomes non-linear the law defined to couple the interface with the CNTs. The formulation is validated studying the elastic response and non-linear behavior of several composites.

In the context of multiscale homogenization, a first-order and an enhanced-first-order formulation is developed. The results obtained for laminate composites using the first-order formulation are compared with other microscopic formulations, showing that the homogenization method is an excellent alternative when microstructural effects must be taken into account. Then, a strategy to conduct non-linear multiscale analysis in an efficient way is proposed. The procedure conserves the dissipated energy through the scales and is mesh independence. The analysis of academic examples is used to show the capacity of the non-linear strategy. Finally, the simulation of an industrial composite component proves the performance and benefits of the non-linear homogenization procedure developed.

# Acknowledgements

This work was financially supported by CIMNE together with the European Community under grant:

NMP-2009-2.5-1 246067 M_RECT “Multiscale Reinforcement of Semi-crystalline Thermoplastic Sheets and Honeycombs”,

FP7-PEOPLE-2013-IRSES 612607 TCAiNMaND “Tri Continental Alliance in Numerical Methods applied to Natural Disasters”, by European Research Council through of Advanced Grant:

ERC-2012-AdG 320815 COMP-DES-MAT “Advanced tools for computational design of engineering materials”, by Dirección General de Investigación Científica y Técnica:

MAT2014-60647-R OMMC “Optimización multi-escala y multi-objetivo de estructuras de laminados compuestos”, by Abengoa Research”, and by Universitat Politecnica de Catalunya (UPC).

All this support is gratefully acknowledged.

# List of Acronyms

 BVP Boundary Value Problem CNTs Carbon Nanotubes CVD Chemical Vapor Deposition DDM Discrete Damage Mechanics EFO Enhanced-First-Order FE Finite Element FEM Finite Element Method FE2 Finite Element Two-Scale FO First-Order IFSS Interfacial Shear Strength LE Linear Element MWCNTs Multiwall Carbon Nanotubes NLAF Non-Linear Activation Function NLS Non-Linear Strategy OpenMP Open Multi-Processing POD Proper Orthogonal Decomposition QE Quadratic Element RHS Right-Hand Side RVE Representative Volume Element SFS Smart First Step SP Serial-Parallel SWCNTs Single Wall Carbon Nanotubes

# Introduction

The continuum mechanics theory has made a great effort to obtain the behavior of homogeneous materials using physical and mathematical concepts showing a good agreement with reality. Furthermore, the constant improvements on computer technology and computer architecture have allowed to improve the numerical tools used to simulate mechanical structures. In the numerical simulation field, one of the most extended methods used for several applications is the Finite Element Method (FEM) [1]. In a FEM analysis the behavior of the homogeneous materials in the structure is simulated by a specific constitutive law or constitutive model with some calibrated parameters.

Composites are non-homogeneous materials formed by two or more different components which can be homogeneous materials or even micro-heterogeneous materials. The homogenized behavior of the composites depends strongly on the internal spatial distribution, the size and the properties of the material components and their respective interfaces. Therefore, composites require more complex and advanced constitutive models than the ones use in single materials.

For a linear analysis or a structural analysis to failure study it is enough to simulate the composite with one orthotropic homogenized characterization and a constitutive law with some complexity. However, more realistic composite constitutive models are necessary to simulate the structures beyond their elastic limit, to obtain the post critic behavior of these or to estimate their tenacity and structural integrity.

The direct application of the FEM is not the most appropriate or effective manner to face the described problem. In a classical FEM analysis each component material has its own constitutive model. Therefore, the numerical model of the structure must to be discretized with a Finite Element (FE) size of at least the size of the components in the composite. In general, this restriction gives as result FE meshes with large number of finite elements which demand an extremely expensive computational cost, and in some cases this analysis is unfeasible to perform. Consequently, to analyze composite structures and to characterize their behavior or fracture modes more suitable strategies must be developed.

## Background and motivation

The complexity of the composite materials has promoted that several formulations appeared to predict their behavior, which are more o less suitable according to the computational cost required, the accuracy in the results desired or even the expected failure type. Further, the development of a new generation of composites with improved properties, more reliable and cheap has extended its use to many industrial applications.

The phenomenological homogenization methods are a possibility to analyze composite materials with a heterogeneous internal structure. In this context, the most usual method is the classical mixing theory proposed initially by Truesdell and Toupin [2]. The formulation obtains the homogenized behavior of the composite through the compatibility equation and from the mechanical performance of the component materials, which are simulated with their own constitutive laws. Later, several modifications and extensions of this classical theory of mixtures have enabled the resolution of any composite with reinforced matrix, without the limitation required by the compatibility equation [3].

One of the most significant modifications of the mixing theory is the Serial-Parallel (SP) continuum approach. In the SP formulation the mechanic characteristic of the composite is obtained using not only the properties and constitutive model of the material components but also taking into account their topological distribution [4]. The SP mixing theory assumes a serial-parallel self-adjusting behavior to the topological distribution of fiber embedded in the matrix of the composite material. This approach imposes the iso-strain condition in the fiber alignment direction on the components of the composite (as parallel materials) and the iso-stress condition in the orthogonal direction (as serial materials).

The reinforcements developed today, as the nanofibers or Carbon Nanotubes (CNTs) and the renewed composites as the reinforced concrete with short fibers, require sophisticated formulations for their simulation. The interface zone between matrix and reinforcement has a meaningful effect in the final properties and response of the composite structures made with these current composites. Many effort have been made to consider the debonding phenomenon in laminated composites but it is not enough to totally characterize the behavior of these micro-heterogeneous composites. Therefore, renovated formulations or renewed modifications of existing theories should be developed to face with the challenge of predict the behavior of these new composites.

The homogenization techniques are another option to analyze composite materials. In these methods the characterization of the entire composite is obtained through the analysis of its internal structure or microscopic structure. In this context, an approach extensively used is the called multiscale homogenization method. In general, the formulation is based on the use of the concept of unit cell or Representative Volume Element (RVE) [5]. The definition of the RVE corresponds with a microscopic subregion which is representative of the entire micro-structure (referred as micro-scale) level of the composite. This is employed to determine the homogenized properties and behavior of the composite level (also known as macro-scale). It is assumed that the RVE must contain a sufficient number of inclusions to make the homogenized moduli independent of homogeneous forces or displacements on the RVE boundary.

Within this context, one of the most extended and popular method is the known first-order homogenization approach [6]. This multiscale method uses the macro-scale deformation gradient tensor (or the strain tensor) to solve the micro-scale problem. The composite behavior (the macro-scale stress-strain relationship) is obtained by a detailed modeling of the internal heterogeneous structure of the composite in the RVE. Therefore, the approach does not require any composite constitutive assumption or compatibility equation to address the composite response. Moreover, there are not restriction about the constitutive law of the component materials, even non-linear materials and time-dependency models can be taken into account. The benefits of the method becomes in a challenge when a non-linear analysis of a three-dimensional structure is studied. Considering a ${\textstyle {\hbox{FE}}^{2}}$ homogenization technique [7], it is required for each time step to solve one RVE at each point of integration at the macro-scale because the non-linear threshold and non-linear behavior of the homogenized composite are unknown. Therefore, the computational cost in the non-linear analysis of an industrial component by using multiscale ${\textstyle {\hbox{FE}}^{2}}$ homogenization is extremely expensive, and in many cases, is unsuitable to perform.

In addition to the computational cost to address the non-linear problem with multiscale homogenization methods, the softening issue must be considered too. The non-linear constitutive law of the component materials are defined in the RVE problem. Consequently, the non-linear behavior starts in the micro-scale and then, it moves up to the macro-scale. Because of this, novel computationally efficient multiscale strategies dealing with non-linear problem should be developed taking into account also the conservation of the dissipated energy through the scales. Besides, they must be macro and micro mesh independence for the case of ${\textstyle {\hbox{FE}}^{2}}$ homogenization.

In the last decade, a second-order computational homogenization was proposed as a natural extension of the first-order homogenization method [8]. It was developed to be applied in critical regions of intense deformation, where the characteristic wave length of the macro-scale deformation field is of the order of the size of the micro-scale. Therefore, in this approach the macroscopic gradient of the deformation gradient is also incorporated in the microscopic scale problem. The first-order equilibrium problem is conserved in the micro-scale though a higher-order equilibrium problem appears in the macro-scale. The solution of the proposed multiscale approach is made through a complex finite element implementation, which restricted its popular application.

The main advantage of the described second-order homogenization is that it can consider intense localization phenomena, then it is a desirable approach for non-linear analysis. On the other hand, the benefit of the first-order homogenization is that it considers first-order equilibrium equations at both scales, which represents an advantage from a point of view of computational implementation. Therefore, an enhanced-first-order approach could be an interesting option to account second-order effects of the macro-scale from the micro-scale by the incorporation of macroscopic second-order deformation measure in the microscopic boundary value problem.

## Objectives

The main objective of this study is to develop a comprehensive formulation for the analysis of three-dimensional composite structures in linear and non-linear range. In this context, the partial targets to address the global aim of this monograph can be written in a synthesized form as:

1. Development of a phenomenological homogenization formulation based on the mixing theory for the analysis of CNTs reinforced composites. The formulation should consider the effect of the CNTs-matrix interface in the composite behavior.
2. Extension to three dimensions of the first-order multiscale homogenization for the numerical analysis of composite structures. Implementation of an elimination of redundant unknowns method to solve the microscopic boundary value problem considering the constraint conditions on the boundary domain.
3. Improvement of the first-order multiscale approach implemented to consider second-order effects in the microscopic scale from the macroscopic scale.
4. Comparison of numerical simulations with other microscopic formulations to show the advantages and drawbacks of the developed first-order multiscale procedure.
5. Development of a non-linear strategy to optimize the computational cost of the analysis of real-size composites structures using a multiscale homogenization approach.

In order to achieve the objectives described previously, there are parallel tasks that must be addressed. These are important milestones of this work that are worth to be mentioned. Among them:

1. Parallel numerical implementation of the developed approaches through an Open Multi-Processing (OpenMP) philosophy in the finite element code PLCd [9].
2. Development of a preprocessor manager to deal with the numerical models of the microstructure (RVE) using a GID problem type [10,11].
3. Numerical validation of the different formulations implemented in PLCd through the simulation of several composites and by the analysis of real-structural components.

## Outline

In the present monograph is possible to observe that from a theoretical point of view the main goal can be divided in two parts. For this reason, the document is arranged in two major self-contained parts.

In Part I a phenomenological homogenization model for the analysis of composites material using CNTs as reinforcement is presented. The formulation developed is based on the mixing theory. In this context, Chapter 2 shows a review of the state of the art of the classical mixing theory and its subsequent modifications while Chapter 3 introduces general considerations about the CNTs and a state of the art of the production methods and of the measured mechanical properties. Then, Chapter 4 presents the formulation and numerical implementation of the “ad hoc” homogenization model developed in this study. Chapter 5 shows the validation and numerical examples analyzed using the formulation proposed in the above chapter. Finally, in Chapter 6 the conclusions and future work about the model developed in this part of the monograph is approached.

In Part 6 the developed multiscale homogenization approach for composite structures is described. The state of the art is addressed in Chapter 8, which shows the fundamental theories and the latest developments about this research topic. Chapter 9 presents the formulations and implementations of the first-order homogenization and the proposed enhanced-first-order extension to consider second-order effects. In Chapter 10 the implemented two-scale homogenization procedure is compared with other micro-structural formulations. Chapter 11 describes the non-linear strategy proposed for multiscale approaches, also its validation and numerical applications are shown. The conclusions and future work of this second part are addressed in Chapter 12.

In the final conclusions chapter the achievements of the present study are exposed along with concluding remarks and future works derived from this monograph.

## Research dissemination

The work included in this monograph resulted in the following scientific publications:

## Part I

1. F. Otero, S. Oller, X. Martinez and O. Salomon. Modelling of behaviour of carbon nanotube-reinforced composites. In: MATCOMP'11 - IX Congreso Nacional de Materiales Compuestos. Girona, España 2011.
2. F. Otero, S. Oller, X. Martinez and O. Salomon. Numerical modelling of behaviour of carbon nanotube-reinforced composites. In: COMPLAS XI - XI International Conference on Computational Plasticity. Fundamentals and Applications. Barcelona, España, 2011.
3. F. Otero, S. Oller, X. Martinez and O. Salomon. Modelling the elastic behavior of carbon nanotube-reinforced composites. In: Composites 2001 - 3rd ECCOMAS Thematic Conference on the Mechanical Response of Composites. Hannover, Germany, 2011.
4. F. Otero, S. Oller, X. Martinez and O. Salomon. Modelling viscoelastic behaviour of carbon nanotube-reinforced thermo-plastics. In: MECOM 2012 - X Congreso Argentino de Mecanica Computacional. Salta, Argentina, 2012.
5. F. Otero, X. Martinez, S. Oller and O. Salomon. Study and prediction of the mechanical performance of a nanotube-reinforced composite. Composite Structures. 2012, 94(9):2920-2930. doi: 10.1016/j.compstruct.2012.04.001.

## Part II

1. X. Martinez, F. Otero and S. Oller. Strategy for an efficient material non-linear multiscale analysis In: COMAT 2015 - VI International Conference on Science and Technology of Composite Materials. Buenos Aires, Argentina, 2015.
2. F. Otero, S. Oller, X. Martinez and O. Salomon. Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Composite Structures. 2015, 122:405-416. doi: 10.1016/j.compstruct.2014.11.041.
3. F. Otero, X. Martinez, S. Oller and O. Salomon. An efficient multi-scale method for non-linear analysis of composite structures. Composite Structures. 2015, 131:707-719. doi: 10.1016/j.compstruct.2015.06.006.
4. F. Otero, S. Oller and X. Martinez. Multiscale computational homogenization: review and proposal of a new enhanced-first-order method. Archives of Computational Methods in Engineering. 2016, 1-27. doi: 10.1007/s11831-016-9205-0

In addition, part of the work was presented at the following unpublished conferences:

1. F. Otero, S. Oller, X. Martinez and O. Salomon. Numerical homogenization for the simulation of composites materials. Comparison with other micro mechanical formulations. Mechanics of Composites (MECHCOMP2014). Long Island, NY State, USA, 8-12 June 2014.
2. F. Otero, S. Oller and X. Martinez. Non-linear multiscale strategy to analyze composite materials efficiently. 18th International Conference on Composites Structures (ICCS18). Lisbon, Portugal 15-18 June, 2015.

Finally, part of the work presented in this document is the result of the collaborating with external researchers during the following research stays:

1. WVU-(West Virginia University), 6-month doctoral research stay. Worked under the direct supervision of Prof. Ever J. Barbero in the Mechanical and Aerospace Engineering of the Benjamin M. Statler College of Engineering and Mineral Resources, WVU. Morgantown, USA. January - June 2014. The following article resulted from the work developed during the stay: M.M. Moure, F. Otero, S.K. García-Castillo, S. Sánchez-Sáez, E. Barbero and E.J. Barbero. Damage evolution in open-hole laminated composite plates subjected to in-plane loads. Composites Structures. 2015, 133:1048-1057. doi: 10.1016/j.compstruct.2015.08.045.
2. CIMAT-(Centro de Investigación en MATemáticas), 2-month research stay in the framework of the TCAiNMaND project, a Marie Curie International Research Staff Exchange Scheme (IRSES) under grant agreement FP7-612607. Worked under the direct supervision of Dr. Salvador Botello in the Computational Sciences Department of CIMAT. Guanajuato, México. November - Dicember 2015.

# 1 Introduction

Composites are materials made of at least two different components. Generally, are constituted by a matrix that surrounds the reinforcing elements, which may be in the form of particles, nanotubes, short fibers, fibers, etc [12]. The main function of matrix component is to give cohesion, support the reinforcement and transfer the external actions to the reinforcements. While the main task of the reinforcement component is to improve the matrix properties. The appropriate design of structural elements made of this type of composite material requires the use of composite constitutive models capable of estimating their stiffness, strength and different failure modes.

In case of using fibers or nanotubes as reinforcement components, the performance of the composite depends on the achievement of the following four main characteristics:

1. Aspect ratio of the fibers. The fiber aspect ratio is a dimensionless geometric measurement that results from dividing the length of the reinforcement by its diameter. This parameter is important because the stress distribution in the reinforcement depends on it [13]. In fibers with high aspect ratios the fiber end effect is less important. The reinforcement is considered to behave as a long fiber when its aspect ratio is larger than 1000.
2. Dispersion of fibers in the composite. A uniform distribution of reinforcement in the composite is fundamental to ensure that it is completely surrounded by matrix. This is necessary to obtain an effective stress transfer. A good dispersion of the reinforcement in the composite also helps to have a more uniform stress distribution in it, reducing the regions susceptible having stress concentrations.
3. Fiber alignment. It has been shown that the difference between random distribution and perfect alignment may represent a factor of five in the composite Young's modulus [14]. Fiber alignment also affects the isotropy of the composite, as perfect alignments increase its anisotropy.
4. Load transmission from the matrix to the fiber. The last and probably most important factor is the interfacial tension between matrix and reinforcement. In general, the loads in a composite structure are introduced through the matrix and are transferred to the reinforcement through the interface [13]. Therefore, the interface can be defined as the region surrounding the reinforcement where this stress transfer takes place. The properties of the composite depend on the properties of this region, and on its ability to transfer the load efficiently.

On the other hand, the external load applied to a composite is shared disproportionately by the different components, as their elastic properties are different. In case of considering an iso-strain hypothesis [2] the stresses on the reinforcement will be larger than in the matrix, as the reinforcement is stiffer than the matrix. This unequal stress distribution generates shear stresses between both materials in a region that can is usually called the interface. The load transfer from the matrix to the reinforcement is produced in this region. Shear stresses in the interface increase proportionally to the external load until a critical value, beyond which the interface breaks. This critical value is known as Interfacial Shear Strength (IFSS) and it limits the stress transfer capacity.

In this context, the classical rule of mixtures was one of the first theories used to address a composite constitutive model from a phenomenological point of view [2]. The theory defines the manner in which all components interact to provide the material performance. The iso-strain hypothesis defined in the mixing theory implies a parallel distribution on the components in the composite. It is possible to think in an inverse mixing theory which replaces the iso-strain assumption by an iso-stress assumption, therefore it means a serial distribution on the components in the composite. The characterization of the composite depends of the hypothesis used in the formulation. Then, modifications of the mixing theory capable to consider all possible behaviors of the composite: parallel, serial and mixed were proposed [15]. Finally, the serial/parallel concept was introduced in the theory, which replaced the iso-strain hypothesis by an iso-strain condition in the fiber direction and iso-stress condition in the transversal directions [4].

The mixing theory began to be considered as a constitutive equation manager when its hypothesis were coupled with a thermodynamical description of the composite components [16]. Therefore, the formulation obtains the relation between the components even when they have reached their elastic limit. With this at hand, the different failure phenomenons present in the composites such as debonding or delamination were modeled through of the constitutive law of the material components.

## 1.1 Part's outline

In this part of the monograph a renewed modification of the mixing theory is proposed to consider the effect of the reinforcement-matrix interface zone in the final response of composite. The present formulation is developed for composites that use CNTs as reinforcement.

Following this aim, Chapter 2 shows an state of the art of the classical mixing theory and its modifications since the early inclusion of the SP concept until the current sophisticated modifications. Chapter 3 presents a review of the state of the art regarding the different production methods and measured mechanical properties found in the literature regarding CNTs is presented. In this chapter is also addressed the different issues that should be considered in a constitutive formulation for reinforced composites with CNTs.

In Chapter 4 the phenomenological homogenization based in the on the mixing theory is developed. The insertion of the concepts of serial and parallel behavior in the CNTs-matrix bonding through of the definition of a parallel factor is shown. The CNTs debonding phenomena is also considered by a material non-linearity proposed. At the end of the chapter, the implementation in the FEM code PLCd is presented.

Chapter 5 shows the results of validation of the implemented composite constitutive model using information data from literature. The validation of the model is made for linear and non-linear behavior using experimental data of several composites. Then, numerical examples are developed showing the good behavior of the phenomenological homogenization model developed.

Finally, Chapter 6 the conclusions about the composite constitutive formulation developed in this part of the study are addressed in detail.

# 2 State of the art

The first part of the present study proposes a phenomenological composite constitutive model which is based on the classical rule of mixtures. Therefore, in the following a literature review is dedicated to explain this mixing theory and the modifications which have been developed over the years.

## 2.1 Classical mixing theory

The mixing theory was originally proposed by Truesdell and Toupin [2]. Later, Truesdell [17] extended the theory to linear systems and Green and Adkins [18] presented a general non-lineal constitutive equations. Finally, Ortiz and Popov [19,20] proposed a general constitutive equations for unreinforced concrete idealized as a composite material.

The classical mixing theory is based on the mechanical of local continuous solid and it is appropriated to explain the behavior of a point in a solid compound. It is based on the principle of interaction of the component substances in the composite material, assuming the following hypothesis: i) on each infinitesimal volume of the composite are involved a set of component substances,ii) each component contributes to the behavior of the compound material in the same proportion as their volumetric fraction, iii) all component materials have the same deformation (compatibility equation or closure equation) and iv) the occupied volume of each component is much smaller than the total volume of the composite.

The second hypothesis implies a homogeneous distribution of all substances in the compound material. The volume fraction, the internal distribution and the interaction between the different substances components, which have its respective constitutive law, determine the behavior of the composite material. This allows combining materials with different behavior (elastic, elasto-plastic, elasto-damage, etc.), which have an evolutionary behavior governed by its own law and internal variables [20,16].

The third hypothesis demands that the following condition of compatibility must be fulfilled

 ${\displaystyle \varepsilon _{ij}=\left(\varepsilon _{ij}\right)_{1}=\left(\varepsilon _{ij}\right)_{2}=\cdots =\left(\varepsilon _{ij}\right)_{n},}$
(2.1)

where the assumption of infinitesimal deformations on each components are considered and where, ${\textstyle \varepsilon _{ij}}$ and ${\textstyle \left(\varepsilon _{ij}\right)_{n}}$ are the strain tensors of the composite and of the ${\textstyle n-th}$ component of the compound material, respectively.

The specific Helmholtz free energy of the composite is given by the sum of the specific Helmholtz free energies of each components of the composite multiplied by its volume fraction, that is

 ${\displaystyle \Psi \left({\boldsymbol {\varepsilon }}^{e},\theta ,{\boldsymbol {\alpha }}\right)=\sum _{c=1}^{n}k_{c}\Psi _{c}\left({\boldsymbol {\varepsilon }}_{c},{\boldsymbol {\varepsilon }}_{c}^{p},\theta ,{\boldsymbol {\alpha }}_{c}\right),}$
(2.2)

where ${\textstyle \Psi _{c}}$ is the specific Helmholtz free energy, ${\textstyle k_{c}}$ is the volume fraction, ${\textstyle {\boldsymbol {\varepsilon }}_{c}^{p}}$ is the plastic strain tensor and ${\textstyle {\boldsymbol {\alpha }}_{c}}$ are the inner variables of each one of the ${\textstyle n-th}$ components in the composite.

The volume fraction coefficient allows to consider the contribution of each material to the composite and it is obtained with the following equation as

 ${\displaystyle k_{c}={\frac {dV_{c}}{dV_{0}}},}$
(2.3)

where ${\textstyle dV_{c}}$ is the volume of the ${\textstyle c-th}$ component and ${\textstyle dV_{0}}$ is the total volume of the composite. The volume fractions of the components must to satisfied the following condition:

 ${\displaystyle \sum _{c=1}^{n}k_{c}=1.}$
(2.4)

Equation 2.4 guarantees the conservation of mass. Following with the procedure used for a simple material [21,22,23,24], from the Clausius-Duhem inequality and applying the Coleman method, the constitutive equation of the composite is obtained as

 ${\displaystyle \sigma _{ij}={\frac {\partial \Psi \left(\varepsilon _{ij},\theta ,\alpha _{i}\right)}{\partial \varepsilon _{ij}}}=\sum _{c=1}^{n}k_{c}{\frac {\partial \Psi _{c}\left(\varepsilon _{ij},\theta ,\alpha _{i}\right)}{\partial \varepsilon _{ij}}}=\sum _{c=1}^{n}k_{c}\left(\sigma _{ij}\right)_{c},}$
(2.5)

where, ${\textstyle \sigma _{ij}}$ and ${\textstyle \left(\sigma _{ij}\right)_{c}}$ are the stress tensors of the composite and of the ${\textstyle c-th}$ component of the compound material, respectively. The composite constitutive tensor is obtained considering the variation of the composite stress tensor respect to the strain tensor, therefore

 ${\displaystyle C_{ijkl}={\frac {\partial \sigma _{ij}}{\partial \varepsilon _{kl}}}={\frac {\partial ^{2}\Psi \left(\varepsilon _{ij},\theta ,\alpha _{i}\right)}{\partial \varepsilon _{ij}\partial \varepsilon _{kl}}}=\sum _{c=1}^{n}k_{c}\left(C_{ijkl}\right)_{c}.}$
(2.6)

However, the closure equation given by 2.1 imposes a strong limitation of the classical theory of mixtures because it is strictly valid only for composites with parallel behavior. Moreover, this limitation is extended to non-linear range because each component can have different deformation for a given load step. Because of this, several formulations have been proposed from the classical mixing theory in order to consider different internal behavior (not only parallel behavior) and the nonlinearity of each components too.

## 2.2 Modifications to the mixing theory

Over the years the classical mixing theory has had many modifications and improvements with the objective of expanding its scope. Some of the most relevant developments are presented in this section.

### 2.2.1 Mixing theory using serial-parallel model

The classical theory of mixtures was modified by Oller et al. [25] and Neamtu et al. [15] introducing the serial-parallel concept. The model allows to represent composites for various possible combinations of serial and/or parallel behavior of their components. The properties of the composite are obtained using the properties of each component and taking into account its topological distribution. The modification is based on the definition of the total strain field as a weighted sum of the contributions of the deformation components in series and parallel. Therefore

 ${\displaystyle \varepsilon _{ij}=\left(1-\aleph \right)\varepsilon _{ij}^{par}+\aleph \,\varepsilon _{ij}^{ser},}$
(2.7)

where ${\textstyle \varepsilon _{ij}}$ is the total strain tensor of the composite, ${\textstyle \varepsilon _{ij}^{par}}$ and ${\textstyle \varepsilon _{ij}^{ser}}$ represent the parallel and serial strain tensor, respectively. And ${\textstyle \aleph }$ is the coupling parameter that relates in weighted form the serial-parallel behavior, it has a value range from 0 to 1. The deformation components in parallel and serial behavior are approximated by

 ${\displaystyle \varepsilon _{ij}^{par}\approxeq {\frac {1}{n}}\sum _{c=1}^{n}\left(\varepsilon _{ij}\right)_{c}\qquad ,\qquad \varepsilon _{ij}^{ser}=\sum _{c=1}^{n}k_{c}\left(\varepsilon _{ij}\right)_{c}.}$
(2.8)

This modification of the classical mixing theory has the disadvantage that the coupling parameter, in general, must be calibrated with experimental tests of the composite.

### 2.2.2 Generalized mixing theory

The proposed modification of the classical theory of mixtures by Oller [3] is a generalization of this theory. The new proposed is enabled to solve any reinforced matrix composite, without the limitation of the compatibility equation (see 2.1) required by the classical theory. The closure equation is satisfied automatically by the proposed modification. The fundamental hypothesis of this theory's generalization is a new definition of the third hypothesis of the classical theory. The new third hypothesis is: iii) the components must satisfy a generalized compatibility equation to fit the topology of the serial-parallel composite. The new hypothesis allows establishing the relationship between the composite deformation and the deformation of each component. The new compatibility equation provides the link between the parallel behavior and the serial behavior and can be expressed as

 ${\displaystyle (\varepsilon _{ij})_{c}=\underbrace {(1-\chi _{c})\cdot I_{ijkl}\,\varepsilon _{kl}} _{\left(\varepsilon _{ij}^{par}\right)_{c}}+\underbrace {\chi _{c}\cdot [(\phi _{ijkl})_{c}\cdot (\varepsilon _{kl}-\varepsilon _{kl}^{p})+(\varepsilon _{kl}^{p})_{c}]} _{\left(\varepsilon _{ij}^{ser}\right)_{c}},}$
(2.9)

where ${\textstyle (\varepsilon _{ij})_{c}}$ is the strain tensor of the ${\textstyle c-th}$ component, which can be separated in its parallel ${\textstyle (\varepsilon _{ij}^{par})_{c}}$ and serial ${\textstyle (\varepsilon _{ij}^{ser})_{c}}$ component, respectively, and ${\textstyle \varepsilon _{ij}}$ is the total strain tensor in the composite. Equation 2.9 can be rewritten as

 ${\displaystyle (\varepsilon _{ij})_{c}=[(1-\chi _{c})\cdot I_{ijkl}+\chi _{c}\cdot (\phi _{ijkl})_{c}]:\varepsilon _{kl}-\chi _{c}({\hat {\varepsilon }}_{kl}^{p})_{c},}$
(2.10)

where ${\textstyle ({\hat {\varepsilon }}_{kl}^{p})_{c}}$ is a plastic strain tensor without physical meaning, which is defined only for operating purpose and it is obtained from the plastic strain tensor of the composite distributed among its components according to ${\textstyle (\phi _{ijkl})_{c}\,\varepsilon _{kl}^{p}}$ and the plastic strain tensor of the current component ${\textstyle (\varepsilon _{kl}^{p})_{c}}$. The serial-parallel coupling parameter is defined as ${\textstyle 0\leq \chi _{c}=\sin \alpha _{\chi }\leq 1}$, where ${\textstyle \alpha _{\chi }}$ corresponds to the existing angle between the reinforcement orientation and the orientation of the higher principal stress.

### 2.2.3 Mixing theory expressed in finite strains

The extension to finite strains of the classical theory of mixtures considers that the third hypothesis (original closure equation given by 2.1) must be verified on the referential configuration and on the spatial configuration for each component[26]

 ${\displaystyle {\begin{array}{l}E_{ij}=(E_{ij})_{1}=(E_{ij})_{2}=\cdots =(E_{ij})_{n}\\\\e_{ij}=(e_{ij})_{1}=(e_{ij})_{2}=\cdots =(e_{ij})_{n},\end{array}}}$
(2.11)

where ${\textstyle E_{ij}}$ is the Green-Lagrange strain tensor and ${\textstyle e_{ij}}$ is the Almansi strain tensor. Considering the definition of the right Gauchy-Green tensor and 2.11 the compatibility equation can be written as a function of the deformation gradient tensor as

 ${\displaystyle F_{ij}=(F_{ij})_{1}=(F_{ij})_{2}=\cdots =(F_{ij})_{n}.}$
(2.12)

The others hypothesis of the classical theory must be also verified. The relationship between the volume of a component in the spatial configuration and in the referential configuration is given by the determinant of the deformation gradient tensor and it is

 ${\displaystyle dV_{c}={\frac {1}{J}}dv_{c}.}$
(2.13)

With 2.13 it is possible to demonstrate that the volume fraction of the components do not change in both configurations.

The solution process starts by estimating the strain increments at the reference configuration and then through tensor transport operations (“push-forward”) the strain tensor in the updated configuration can be obtained. The constitutive equation of each components of the composite is integrated in the updated configuration. Each of these components may have different kinds of constitutive behavior (plasticity, damage, etc.) and also, these constitutive models may be isotropic or anisotropic. Then, with the integrated stress state of the components it is possible to obtain the stress state and the constitutive tensor of the composite. Finally, the obtained composite informations are transported (“back-forward”) to the reference configuration and then, the internal forces are computed. The balance between the internal forces and applied external forces is verified in an iterative procedure until convergence.

### 2.2.4 Generalized mixing theory expressed in finite strains

The procedure to extend the generalized theory of mixtures to finite strains is the same than the one used to extend the classical mixing theory in Section 2.2.3. This generalized theory starts with the non-compliance of the classical compatibility equation. Therefore, the proposed new closure equation given by 2.10 must be written now in the reference and the updated configuration, that is

 ${\displaystyle {\begin{array}{c}(E_{ij})_{c}=[(1-\chi _{c})\cdot I_{ijkl}+\chi _{c}\cdot (\Phi _{ijkl})_{c}]:E_{kl}-\chi _{c}({\hat {E}}_{kl}^{p})_{c}\\\\(e_{ij})_{c}=[(1-\chi _{c})\cdot I_{ijkl}+\chi _{c}\cdot (\phi _{ijkl})_{c}]:e_{kl}-\chi _{c}({\hat {e}}_{kl}^{p})_{c}.\end{array}}}$
(2.14)

Equations 2.14 give the strain tensors for each component in both configurations. The constitutive equations are obtained following a similar formulation used in the classical theory of mixtures in finite strains. Finally, the stress tensors of the composite in the reference configuration ${\textstyle {S}}$ and in the current configuration, the Kirchoff stress tensor ${\textstyle {\boldsymbol {\tau }}}$ are, respectively

 ${\displaystyle \mathbf {S} =\sum _{c=1}^{n}k_{c}{\scriptstyle [(1-\chi _{c})\cdot {I}_{4}+\chi _{c}\cdot (\Phi )_{c}]^{T}:[({C}^{S})_{c}:[(1-\chi _{c})\cdot {I}_{4}+\chi _{c}\cdot (\Phi )_{c}]:({E}^{e})_{c}]}}$
(2.15)
 ${\displaystyle {\boldsymbol {\tau }}=\sum _{c=1}^{n}k_{c}{\scriptstyle [(1-\chi _{c})\cdot {I}_{4}+\chi _{c}\cdot (\Phi )_{c}]^{T}:[({c}^{\tau })_{c}:[(1-\chi _{c})\cdot {I}_{4}+\chi _{c}\cdot (\Phi )_{c}]:({e}^{e})_{c}]}=J{\boldsymbol {\sigma }},}$
(2.16)

where ${\textstyle {I}_{4}}$ is the fourth order identify tensor, ${\textstyle ({C}^{S})_{c}}$ and ${\textstyle ({c}^{\tau })_{c}}$ are the tangent constitutive tensors for the ${\textstyle c-th}$ component in each configuration, ${\textstyle ({E}^{e})_{c}}$ and ${\textstyle ({e}^{e})_{c}}$ are the elastic strain tensors and ${\textstyle {\boldsymbol {\sigma }}}$ is the Cauchy stress tensor.

### 2.2.5 Mixing Theory by short fiber reinforcements

The formulation of the mixing theory is oriented to a composite where the reinforcements are long fibers, and the condition of the compatibility equations is verified. However, when the aspect ratio of the fiber decreases, the condition of fiber-matrix compatibility is not satisfied. This is because the effect of slip and the limit transmission of forces between fiber and matrix at the ends of the fiber take increasingly significant. This situation creates conditions of stress concentration and distortion in the fiber and the surrounding matrix because of the discontinuity. The effectiveness of the fibers in the composite stiffness decreases when the length of the fiber decreases.

Equation 2.17 shows the axial stress distribution along the fiber [13] as

 ${\displaystyle \sigma _{f}(x)=C_{f}^{\sigma }E_{m}\left[1-{\frac {\cosh \left(\beta \left({\frac {l}{2}}-x\right)\right)}{\cosh \left(\beta {\frac {l}{2}}\right)}}\right]\qquad \qquad \forall \quad 0\leq x\leq {\frac {l}{2}},}$
(2.17)

where ${\textstyle C_{f}^{\sigma }}$ is the Young's modulus, ${\textstyle l}$ is the length of the reinforcement, ${\textstyle E_{m}}$ is the longitudinal strain of the matrix and the parameter ${\textstyle \beta }$ is defined as

 ${\displaystyle \beta ={\sqrt {{\frac {G_{c}}{C_{f}^{\sigma }}}{\frac {2\pi }{A_{f}ln{\frac {r^{'}}{r}}}}}},}$
(2.18)

where ${\textstyle A_{f}}$ is the cross section of the fiber, ${\textstyle G_{c}}$ is the shear modulus of the composite and ${\textstyle r^{'}}$ is the mean distance between the reinforcing fibers.

One way to consider the contribution of the short fiber reinforcement in the classical mixing theory is through the average stress along the fiber, then

 ${\displaystyle {\bar {\sigma }}_{f}={\frac {1}{l}}\int _{0}^{l}\sigma _{f}(x)\,dx=C_{f}^{\sigma }\left[1-{\frac {\tanh \left(\beta {\frac {l}{2}}\right)}{\left(\beta {\frac {l}{2}}\right)}}\right]E_{m}={\widetilde {C}}_{f}^{\sigma }E_{m}.}$
(2.19)

Here, ${\textstyle {\widetilde {C}}_{f}^{\sigma }}$ is the average or homogenized Young's modulus of the reinforcement, which is function of the length of the fiber and of the geometric parameters of the composite.

The obtained short fiber homogenized Young's modulus is smaller than the real fiber Young's modulus, this shows that its participation on the mechanical properties of the composite depend not only of its mechanical properties but also of the overall properties of the matrix-reinforcement assemblage. The same concept used to homogenize the stress along the fiber can be extended to get the three dimension homogenized constitutive tensor of the short fiber reinforcement as

 ${\displaystyle {\widetilde {C}}_{f}^{S}=C_{f}^{S}\left[1-{\frac {\tanh \left(\beta {\frac {l}{2}}\right)}{\left(\beta {\frac {l}{2}}\right)}}\right],}$
(2.20)

where ${\textstyle C_{f}^{S}}$ is the orthotropic constitutive tensor in the referential configuration of the reinforcement. Using the previously described concept, the incorporation of the short fiber in the theory of mixtures can be extend to finite strains too [3].

### 2.2.6 Serial-Parallel (SP) continuum approach

The SP continuum approach has been proposed by Rastellini et al. [4], and it is a natural evolution of the parallel mixing theory developed by Car et al. [26,3]. The theory is based on the compatibility conditions defined by Trusdell and Toupin [2], but introduces a modification in the iso-strain hypothesis. The iso-strain condition is imposed in the reinforcement direction (normally fiber) and a new iso-stress condition is imposed in the transversal directions. The theory is based on the following hypotheses:

• The constituent materials of the composite are subjected to the same strain in the parallel (fiber) direction.
• Constituent materials are subjected to the same stress in the serial direction.
• The response of the composite material is directly related to the volume fractions of its constituent materials.
• The phases in the composite are considered to be homogeneously ditributed.
• The constituent materials are considered to be perfectly bonded.

In this formulation the definition of the constitutive model for the composite needs the introduction of additional equations that specify somehow the interaction between the component phases. Then, the resulting composite material model will depend crucially on the adopted specific additional equations that characterizes the mechanical interaction at the micro-scale. These additional sets of equations are referred to as “closure equations” and are obtained taking the iso-strain hypothesis in the reinforcement direction and iso-stress hypothesis in transversal directions. Considering only two composite components, the equations that define the stress (${\textstyle \sigma }$) equilibrium and setting up the strain (${\textstyle \varepsilon }$) compatibility between the individual components follow the hypothesis previously described are:

• Parallel behavior
 ${\displaystyle {\begin{array}{l}^{c}\varepsilon _{p}={^{m}\varepsilon }_{p}={^{f}\varepsilon }_{p}\\^{c}\sigma _{p}={^{m}k}^{m}\sigma _{p}+{^{f}k}{^{f}\sigma _{p}}\end{array}}}$
(2.21)
• Serial behavior
 ${\displaystyle {\begin{array}{l}^{c}\varepsilon _{s}={^{m}k}^{m}\varepsilon _{s}+{^{f}k}{^{f}\varepsilon }_{s}\\^{c}\sigma _{s}={^{m}\sigma }_{s}={^{f}\sigma }_{s}\end{array}}}$
(2.22)

where, the superscripts ${\textstyle c}$, ${\textstyle m}$ and ${\textstyle f}$ stand for composite, matrix and fiber, respectively, the subscripts ${\textstyle s}$ and ${\textstyle p}$ correspond to the serial and parallel behavior and ${\textstyle ^{i}k}$ is the volume fraction of each constituent in the composite.

Composite materials that can be modeled with this formulation are those formed of long fibers embedded in a matrix. The theory predicts the different behavior of the composite, depending on the load direction. This formulation can obtain the linear and non linear behavior of structural elements made of composite materials as has been proved in several papers [27,28,29,30,31,32]. The SP theory is able to simulate the delamination problem naturally, without having to define specific elements or predefine the path of fracture. The approach has been also extended to tri-dimensional framework by Martinez et al. [33] and applied for the numerical simulation of structures of reinforced concrete retrofitted with carbon fiber reinforced polymers. In this study the tangent constitutive tensor of each component of the composite is obtained by means of a perturbation method.

# 3 Carbon nanotubes reinforced polymers

Since their discovery by Lijima in 1991 [34], CNTs are considered a new generation of reinforcements [35]. Their “nano” size structure makes them potentially free of defects, which gives them excellent physical properties [36,37]. A nanotube is a tubular cylinder formed by sp2 bonds between the carbon atoms along its length. There are two main nanotube types: Single Wall Carbon Nanotubes (SWCNTs), which are made of a single wall tube with an outer diameter in the order of 1 nm; and Multiwall Carbon Nanotubes (MWCNTs), which consist in several concentric walls, one inside the other, separated by a distance of 0.34 nm [34]. The diameters range of MWCNT varies from 2 to 100 nm. MWCNT can have lengths up to 100 ${\textstyle \mu }$m.

Carbon nanotubes can be obtained by several procedures. The first method used was the arc-discharge [38], which consists in generating an arc discharge between two graphite electrodes in an inert gas atmosphere at low pressure. The continuous electric discharge sublimates the carbon atoms of the electrodes and forms a plasma around them. This method produces free defect nanotubes along their length. The length of these nanotubes can reach 50 ${\textstyle \mu }$m. Another procedure is the laser ablation. This consists in vaporizing the graphite by radiation with a laser pulse, in an inert gas atmosphere, inside a high temperature reactor. The nanotubes are formed when the graphite vapor touches the cold walls of the reactor. Finally, the most common procedure used for commercial production of carbon nanotubes is the deposition of Catalytic Vapour Phase (also named, Chemical Vapor Deposition (CVD)). This procedure allows producing large amounts of nanotubes at a low cost. This method prepares a substrate with a metal layer. The nanotube diameter depends on the size of the metal particles. The process starts by mixing two gases; one of them is used as a source of carbon, and the other for the process itself. The nanotubes grow on the side of the metal catalyst. The generated nanotubes have defects on its surface. This method can provide oriented nanotubes if there is plasma during their growth.

Nanotubes obtained by arc-discharge have Young's modulus values in the order of 1TPa. Recent measurements carried out in arc-MWCNTs (multiwall nanotubes made by arc-discharge) have provide Young's modulus values with values varying from 0.27 to 0.95 TPa, ultimate strain values higher than 12%, and ultimate tensile stresses in the range of 11 to 63 GPa [39]. In these measurements it was also obtained the stress-strain curve of the MWCNTs with help an electric microscope.

The properties obtained for CVD-MWCNTs (multiwall carbon nanotubes obtained by CVD) are low due to the defects in the nanotubes surface. The firsts Young modulus measurement known was made with an atomic force microscope [40] and the values obtained were in the range of 12 to 50 GPa. Later on, new measurements have shown Young modulus values in order of 0.45 TPA, and ultimate tensile stresses of 3.6 GPa [41]. The lower measured values were associated with defects in the nanotube and with the slipping of the inner tubes in MWCNTs. The difference in measured values between CVD-MWCNTs and arc-MWCNTs shows the influence of defects on the properties of these new materials.

It is not entirely clear which nanotube type performs better as a reinforcement. A recent study made by Cadeck et al. [42] comparing the properties of a polyvinylalcohol (PVA) matrix reinforced with different types of CNTs nanotubes (double wall nanotubes (DWCNT), SWCNTs, arc-MWCNTs and CVD-MWCNTs) showed that the effectiveness of reinforcement is inversely proportional to its diameter, except when using SWCNTs. The study also proved that the composite properties are proportional to the total interface area. The composite reinforced with SWCNTs had the lowest properties; this result is associated with slipping of SWCNTs inside the bundles. Finally, the study states that the best properties are obtained with the CVD-MWCNTs with smaller diameter.

Currently, there are several methods that can be used to produce nanotube-reinforced composites. The choice of the most appropriate method depends of nature of the involved components [35]. All methods seek to produce a composite with a good dispersion of the CNT reinforcement and to create an interface capable of transmitting the external load to the nanotubes. The manufacturing process has to be selected taking into account that it must not affect the properties of the composite components.

Several studies have shown that the composite formation generates an interface zone around the carbon nanotubes. This interface has a different morphology and properties than the original matrix [43,44]. The size, shape and properties of the interface have a strong dependence on the matrix type [45] and the formation process. Pull out experimental tests indicate that IFSS values are much higher than the theoretical ones [46], which are calculated using the shear strength of the matrix. This result suggests that the interface region around the nanotube has better properties than the rest of the matrix [47]. Some studies estimate that in this region matrix properties may improve by an order of magnitude [48]. Fracture surface images obtained from composites with strong nanotube-matrix bond show that the interface zone has a thickness several times larger than the nanotube diameter [49,50]. In the case of semi-crystalline matrices, the interface zone is associated with crystal nucleation around the nanotubes [51].

All manufacture processes seek to obtain a composite with a strong bond between the nanotube and the matrix, in order to transfer effectively the loads. The IFSS defines the capacity of the bond. Values of 500 MPa have been obtained for the IFSS when observing the stresses induced to a broken nanotube, these values where obtained using a Transmission Electron Microscope (TEM). The study attributes this value to the presence of covalent bonds between the matrix and the nanotube [52]. Molecular Dynamics (MD) simulations carried out confirm that strong bonds are obtained when these are covalent. In fact, the transfer load of the interface increases by an order of magnitude with just a 1% of covalent bonds in its surface [53]. On the other hand, the generation of many covalent bonds in the interface is detrimental to the intrinsic properties of the nanotube [53,54].

When there are not covalent bonds, the interaction between matrix and nanotube is made with Van der Waals forces. Several studies show that this union is weaker. Molecular Dynamics simulations made by [53] predicted values of the IFSS that do not exceed 2.8 MPa. Another study made by [55] predicted values up to 160 MPa. According to [56], the differences in the results depend on the polymer type and they can be in the range of 80 to 135 Mpa. The difference in the results, and the good values of IFSS, were attributed to the morphology and the capacity of the matrix to generate helical chains around the nanotube. On the other hand, nanotubes have a smoother outer surface and therefore, the contribution of the frictional forces to the IFSS are an order of magnitude lower [57].

Experimental results of pull-out tests show values of IFSS between 20-90 MPa [48,57]. Other experiments using the drag-out technique have shown values between 35-376 MPa [47]. The disparity of the results suggests that is not always possible to generate covalent bonds. The maximum values obtained experimentally are associated to covalent bonds and consider that the interface zone has better properties than the rest of the matrix.

Carbon nanotubes, mainly SWCNTs, tend to agglomerate. This makes very difficult to obtain a good dispersion of those in the polymer. Besides, the smooth surface of the nanotubes leads to a possible lack of bond between the nanotube and the matrix. Currently these problems are solved with a chemical functionalization of the CNTs. The covalent functionalization can be done by modifying the carboxylic acid groups on the nanotube surface and or by direct addition of reagents. The drawback of functionalizing the nanotubes is that there is an intrinsic degradation of their properties [54]. In general, two different methods have been used for the functionalization: “grafting from” and “grafting to”.

The “grafting from” method is based on the initial immobilization of initiators on the nanotube surface, followed by an in situ polymerization of the suitable matrix for the formation of polymer molecules around the nanotube [58,59]. The advantage of this method is that it allows the formation of composites with a high density of nanotubes. The disadvantage is that this method requires strict control of the quantities and the conditions in which the polymerization reaction takes place.

The “grafting to” method makes the union of preformed polymer molecules to functional groups on the surface of the nanotube through chemical reactions [60,61]. The advantage of this method is that it can be used with commercial polymers. However, it has as a limitation that the initial union of the polymer chains inhibits the diffusion of macromolecules to the surface. Therefore, the density of functionalization is low.

The above description shows that the final properties of the composite depend on many parameters. Together with these, there are others aspects that may also condition the final properties of the composite, such as the undulation and misalignment of the nanotubes inside the matrix. All this variability can be considered the responsible of not having yet an accepted theory capable of describing correctly the performance of nanotube-reinforced composites. It is also the reason because the existing theories fail in their predictions. Comparisons between measured mechanical properties and theoretical results, show that the theoretical predictions are generally three times higher than measured results [50,62].

# 4 Phenomenological homogenization of CNTs reinforced polymers

Carbon nanotubes have been regarded as ideal reinforcements of high performance composites. A key factor for the reinforcement efficiency is the interface bonding between the CNTs and the matrix. In this chapter the formulation and numerical implementation of a new constitutive model to predict the performance of composites made of CNTs is presented [63]. The composite constitutive model takes intro account explicitly the mechanical performance of the interface between the matrix and the CNTs. The proposed model is based in the classical mixing theory. As it is written, the mixing theory can be understood as a constitutive model manager. Therefore, the mechanical performances of the composite are obtained from the behavior of the composite components, each one simulated with its own constitutive law [26,64]. The present new composite constitutive model is formulated with the same philosophy, which increases its versatility and simulation capability.

## 4.1 Description of the composite constitutive model

The proposed composite constitutive model assumes that the composite is a combination of three different materials: matrix, CNTs and an interface [65]. The interface component corresponds to the matrix that surrounds the CNTs. It is considered as an independent component, with its own constitutive law. The interface is used to define the capacity of the matrix to transfer the loads to the reinforcement.

Although the phenomenological performance of the composite already justifies the definition of an interface material; images obtained with Scanning Electron Microscope (SEM) of CNTs reinforced composites, such the ones shown in Figure 1, prove its actual existence. These images reveal that the structures protruding from the fractured surface have larger diameters than the original MWCNTs used in the sample preparation [49]. The material surrounding the CNTs corresponds to the interface. The presence of an interface, as a differentiable material, is also proved by Differential Scanning Calorimetry (DSC) measurements carried out in composites with a semi-crystalline polymer as matrix. These measurements show a linear increase of crystalline matrix as the nanotube volume fraction increases, suggesting that each nanotube has a crystalline coating [66].

 Figure 1: SEM image of nanomanipulation and fracture surface of composites [49].

Once having conveyed the necessity of including the interface material in the formulation to simulate the mechanical performance of CNTs reinforced composites, in the following is described the new procedure proposed, which is summarized in Figure 2. This figure shows that the composite is divided in several layers, each one containing carbon nanotubes with a different orientation. All layers are coupled together using the parallel mixing theory. This is, assuming that all layers have the same deformation. The new formulation developed provides the mechanical performance of each layer by combining the response of the three coexisting materials: matrix, interface and CNTs. The layer response depends on the materials and on their volumetric participation in the composite.

 Figure 2: Representation of formation for reinforced composite.

First, the layer is split into matrix and a new material that results of coupling the CNTs with the interface. The relation between the matrix and the CNT-interface material is established in terms of the parallel mixing theory (they are assumed to have an iso-strain behavior). On the other hand, CNTs and the interface are coupled together with a combination of parallel and serial mixing theories. The serial mixing theory assumes that all components have the same stresses.

Figure 3 shows scheme used to obtain the performance of the CNT-interface material. This is based in the short-fiber model developed by Jayatilaka [13]. According to this model, the load is transferred from the interface to the nanotube at the ends of the reinforcement, through shear stresses. In this region normal stresses in the fiber increase from zero to their maximum value, which is reached in the central part of the reinforcement. In this region there is not load transfer and shear stresses are null. This whole stress transfer scheme can be simplified assuming a CNT-interface performance defined by a serial mixing theory at the ends of the reinforcement and a parallel mixing theory at the center of it.

 Figure 3: Different regions in the new material CNT-interface.

A parallel factor named ${\textstyle N^{par}}$ is defined to differentiate these two regions. This parameter, multiplied by the nanotube length, provides the length of the nanotube-interface element with a parallel behavior. The length with a serial performance is defined by the complementary factor.

## 4.2 Formulation of the composite constitutive model

The Helmholtz free energy [67] of a material point subjected to infinitesimal deformations can be described with the following thermodynamic formulation [16,22]

 ${\displaystyle \Psi =\Psi \left({\boldsymbol {\varepsilon }},\theta ,{\boldsymbol {\alpha }}\right),}$
(4.1)

where ${\textstyle {\boldsymbol {\varepsilon }}}$ is the strain tensor, ${\textstyle \theta }$ a measure of temperature and ${\textstyle {\boldsymbol {\alpha }}=\left\{{\boldsymbol {\varepsilon }}^{p},d,s\right\}}$ a set of inner variables, for example: ${\textstyle {\boldsymbol {\varepsilon }}^{p}}$ is the plastic strain tensor, ${\textstyle d}$ damage inner variable and ${\textstyle s}$ any other material internal variables.

The proposed model simulates the composite combining the different components using the serial and parallel mixing theories. If this combination is performed according to what has been described in previous Section 4.1, the expression of the Helmholtz free energy may be written as

 ${\displaystyle {\begin{array}{rl}\Psi &={\textstyle k_{m}\Psi _{m}}+\\&\\&+{\textstyle \left(k_{nt}+k_{iz}\right)}[{\underset {{\tilde {\Psi }}_{ntiz}^{par}}{\underbrace {N^{par}{\scriptstyle \left({\overline {k}}_{nt}\Psi _{nt}+{\overline {k}}_{iz}\Psi _{iz}\right)}} }}+{\underset {{\tilde {\Psi }}_{ntiz}^{ser}}{\underbrace {\left(1-N^{par}\right){\scriptstyle \left({\overline {k}}_{nt}\Psi _{nt}+{\overline {k}}_{iz}\Psi _{iz}\right)}} }}],\end{array}}}$
(4.2)

where ${\textstyle \Psi _{m}}$, ${\textstyle \Psi _{nt}}$ and ${\textstyle \Psi _{iz}}$ are the specific Helmholtz free energy for the matrix, the nanotube and the interface components, respectively; ${\textstyle k_{m}}$, ${\textstyle k_{nt}}$ and ${\textstyle k_{iz}}$ are the volume fraction of each component, ${\textstyle N^{par}}$ is the parallel factor and,

 ${\displaystyle {\overline {k}}_{nt}={\frac {k_{nt}}{k_{nt}+k_{iz}}}\qquad {\overline {k}}_{iz}={\frac {k_{iz}}{k_{nt}+k_{iz}}}}$
(4.3)

are the volume fractions of the carbon nanotubes and the interface in the new CNT-interface material. These volume fractions must verify

 ${\displaystyle k_{m}+k_{nt}+k_{iz}=1\qquad {\overline {k}}_{nt}+{\overline {k}}_{iz}=1.}$
(4.4)

The relation among the strain tensors of the different components is

 ${\displaystyle {\boldsymbol {\varepsilon }}={\boldsymbol {\varepsilon }}_{m}={\boldsymbol {\varepsilon }}_{ntiz}^{par}={\boldsymbol {\varepsilon }}_{ntiz}^{ser},}$
(4.5)

being ${\textstyle {\boldsymbol {\varepsilon }}}$ and ${\textstyle {\boldsymbol {\varepsilon }}_{m}}$ the composite and matrix strain tensor, respectively; ${\textstyle {\boldsymbol {\varepsilon }}_{ntiz}^{par}}$ the strain tensor of the new CNT-interface material with a parallel behavior; and ${\textstyle {\boldsymbol {\varepsilon }}_{ntiz}^{ser}}$ the strain tensor of the CNT-interface material with a serial behavior.

The tangent constitutive tensor of the composite material may be derived from 4.2 as

 ${\displaystyle \mathbf {C} ={\frac {\partial ^{2}\Psi }{\partial {\boldsymbol {\varepsilon }}\otimes \partial {\boldsymbol {\varepsilon }}}}=k_{m}{\frac {\partial ^{2}\Psi _{m}}{\partial {\boldsymbol {\varepsilon }}_{m}\otimes \partial {\boldsymbol {\varepsilon }}_{m}}}+{\frac {\partial ^{2}{\tilde {\Psi }}_{ntiz}^{par}}{\partial {\boldsymbol {\varepsilon }}_{ntiz}^{par}\otimes \partial {\boldsymbol {\varepsilon }}_{ntiz}^{par}}}+{\frac {\partial ^{2}{\tilde {\Psi }}_{ntiz}^{ser}}{\partial {\boldsymbol {\varepsilon }}_{ntiz}^{ser}\otimes \partial {\boldsymbol {\varepsilon }}_{ntiz}^{ser}}}.}$
(4.6)

A parallel behavior means that all composite constituents have the same strain value. Therefore:

 ${\displaystyle {\boldsymbol {\varepsilon }}_{ntiz}^{par}={\boldsymbol {\varepsilon }}_{nt}={\boldsymbol {\varepsilon }}_{iz}}$
(4.7)

 ${\displaystyle {\frac {\partial ^{2}{\tilde {\Psi }}_{ntiz}^{par}}{\partial {\boldsymbol {\varepsilon }}_{ntiz}^{par}\otimes \partial {\boldsymbol {\varepsilon }}_{ntiz}^{par}}}=N^{par}\left[{\overline {k}}_{nt}\mathbf {C} _{nt}+{\overline {k}}_{iz}\mathbf {C} _{iz}\right]=N^{par}\mathbf {C} _{ntiz}^{par}.}$
(4.8)

And, a serial behavior means that all composite constituents have the same stress value. Thus:

 ${\displaystyle {\begin{array}{c}{\boldsymbol {\sigma }}_{ntiz}^{ser}={\boldsymbol {\sigma }}_{nt}={\boldsymbol {\sigma }}_{iz}\\\\{\boldsymbol {\varepsilon }}_{nt}=\mathbf {C} _{nt}^{-1}:\mathbf {C} _{ntiz}^{ser}:{\boldsymbol {\varepsilon }}_{ntiz}^{ser}\quad ;\quad {\boldsymbol {\varepsilon }}_{iz}=\mathbf {C} _{iz}^{-1}:\mathbf {C} _{ntiz}^{ser}:{\boldsymbol {\varepsilon }}_{ntiz}^{ser}\end{array}}}$
(4.9)

 ${\displaystyle {\frac {\partial ^{2}{\tilde {\Psi }}_{ntiz}^{ser}}{\partial {\boldsymbol {\varepsilon }}_{ntiz}^{ser}\otimes \partial {\boldsymbol {\varepsilon }}_{ntiz}^{ser}}}=\left(1-N^{par}\right)\left[{\overline {k}}_{nt}\mathbf {C} _{nt}^{-1}+{\overline {k}}_{iz}\mathbf {C} _{iz}^{-1}\right]^{-1}=\left(1-N^{par}\right)\mathbf {C} _{ntiz}^{ser}.}$
(4.10)

Replacing 4.8 and 4.10 in 4.6 it is possible to obtain a simplified expression of the tangent constitutive tensor as

 ${\displaystyle \mathbf {C} =k_{m}\mathbf {C} _{m}+\left(k_{nt}+k_{iz}\right)\left[N^{par}\mathbf {C} _{ntiz}^{par}+\left(1-N^{par}\right)\mathbf {C} _{ntiz}^{ser}\right].}$
(4.11)

The formulation developed require all composite components to fulfill the expression given by 4.1. Therefore, it is possible to use any constitutive law to describe the mechanical performance of the different components.

### 4.2.1 Definition of the parallel factor

The parallel factor is defined as

 ${\displaystyle N^{par}={\frac {l_{par}}{l_{nt}}}\qquad ,\qquad {0}\leq N^{par}\leq {1},}$
(4.12)

where ${\textstyle l_{nt}}$ is the length of the nanotube and ${\textstyle l_{par}}$ is function of geometry and mechanical properties of the nanotube and the interface. The value of this length can be obtained from the equation of tension distribution in a reinforcement considering perfect bond with the matrix, which is [13]

 ${\displaystyle {\boldsymbol {\sigma }}_{nt}\left(x\right)=E_{nt}\left[1-{\frac {\cosh \left(\beta \left(l_{nt}-2x\right)\right)}{\cosh \left(\beta l_{nt}\right)}}\right]{\boldsymbol {\varepsilon }}_{m}}$
(4.13)
 ${\displaystyle \beta ={\sqrt {\frac {2G_{iz}}{E_{nt}d_{nt}^{2}\ln \left(1+{\frac {b}{r_{nt}}}\right)}}},}$
(4.14)

where ${\textstyle x}$ represents the longitudinal positions in the reinforcement, and the subscripts “nt” and “iz” refers to the properties of nanotube and interface zone, respectively. ${\textstyle E}$ and ${\textstyle G}$ are the Young's modulus and the shear modulus, and ${\textstyle b}$ is the thickness material around of the CNTs associated with the interface zone.

Defining ${\textstyle l_{par}=l_{nt}-2x}$, its value can be obtained by finding the position “x” for which the effective modulus obtained from the integration of the tension distribution becomes

 ${\displaystyle E_{eff}={\frac {l_{par}}{l_{nt}}}E_{ntiz}^{par}+\left(1-{\frac {l_{par}}{l_{nt}}}\right)E_{ntiz}^{ser}.}$
(4.15)

This procedure provides a value of the parallel length of

 ${\displaystyle l_{par}={\frac {1}{\beta }}\cosh ^{-1}\left[{\frac {1}{3}}\cosh \left(\beta l_{nt}\right)\right].}$
(4.16)

### 4.2.2 Definition of the volume fraction of the interface region

Based on the results reported in [66], the interface zone can be considered the region surrounding the carbon nanotube in which an amorphous matrix becomes crystalline. The volume fraction of the interface zone can be obtained as

 ${\displaystyle \chi _{c}=\chi _{o}+k_{iz},}$
(4.17)

where ${\textstyle \chi _{c},\,\chi _{o}}$ are the volume fractions of crystalline matrix with and without CNTs, respectively. Assuming that the interface zone is a cylinder around the CNTs, it is possible to relate the volume fraction of the interface zone with the parameter ${\textstyle {\frac {b}{r_{nt}}}}$ as

 ${\displaystyle {\textstyle k_{iz}={\frac {N\left(\pi r^{2}l_{nt}-\pi r_{nt}^{2}l_{nt}\right)}{V}}={\frac {N\pi r_{nt}^{2}l_{nt}}{V}}\left[\left({\frac {r}{r_{nt}}}\right)^{2}-1\right]=k_{nt}\left[\left({\frac {r}{r_{nt}}}\right)^{2}-1\right]},}$
(4.18)

where ${\textstyle V}$ is the total composite volume, ${\textstyle r}$ is the radius of interface zone and ${\textstyle N}$ is the total number of nanotubes in the composite.

The relation between the radius of the nanotube and the interface is obtained replacing 4.18 in 4.17 as

 ${\displaystyle {\frac {r}{r_{nt}}}=+{\sqrt {{\frac {\left(\chi _{c}-\chi _{o}\right)}{k_{nt}}}+1}}\qquad ,\qquad \chi _{c}\geq \chi _{o},}$
(4.19)

and therefore

 ${\displaystyle {\frac {r}{r_{nt}}}=1+{\frac {b}{r_{nt}}}\quad \Rightarrow \quad {\frac {b}{r_{nt}}}=+{\sqrt {{\frac {k_{iz}}{k_{nt}}}+1}}-1\qquad ,\qquad k_{iz}\geq {0.}}$
(4.20)

### 4.2.3 Equivalent properties for MWCNTs

MWCNTs consist of concentric SWCNTs joined together with relatively weak van der Waals forces. For this reason, the capacity to transfer the load from the external wall to the internal walls is low. Some papers [68,69] propose to simulate the CNTs like a solid cylinder with same exterior diameter and length, but with effective properties. The effective properties are obtained assuming that the outer wall takes the total load. In this approach it is assumed that the properties of the outer wall correspond to those of a graphite sheet. The effective stiffness of the MWCNT is calculated by imposing that for a same applied force, the deformation must be the same

 ${\displaystyle {\bar {\boldsymbol {\varepsilon }}}_{nt}={\boldsymbol {\varepsilon }}_{nt}\quad \Rightarrow \quad {\bar {E}}_{nt}={\frac {A_{ow}}{{\bar {A}}_{nt}}}E_{g},}$
(4.21)

where ${\textstyle {\bar {E}}_{nt}}$ and ${\textstyle E_{g}}$ are the Young's modulus of the effective solid nanotube and graphite sheet, respectively, and ${\textstyle {\bar {A}}_{nt}}$ and ${\textstyle A_{ow}}$ are the areas of the effective solid nanotube and outer wall, respectively. Equation 4.21 can be also read as

 ${\displaystyle {\bar {E}}_{nt}=\left[1-\left(1-{\frac {2t}{d_{nt}}}\right)^{2}\right]E_{g}\qquad ,\qquad {\frac {t}{d_{nt}}}\leq {0.5},}$
(4.22)

being ${\textstyle t}$ the thickness of one wall in the MWCNT and ${\textstyle d_{nt}}$ is the external diameter of the MWCNT.

Using the same procedure it is possible to obtain the shear modulus of the solid cylinder, by forcing the same twist when applying the same torque (T).

 ${\displaystyle {\bar {\phi }}_{nt}=\phi _{nt}\quad \Rightarrow \quad {\frac {Tl_{nt}}{{\bar {G}}_{nt}{\bar {J}}_{nt}}}={\frac {Tl_{nt}}{G_{g}J_{ow}}}\quad \Rightarrow \quad {\bar {G}}_{nt}={\frac {J_{ow}}{{\bar {J}}_{nt}}}G_{g},}$
(4.23)

where ${\textstyle {\bar {G}}_{nt}}$ and ${\textstyle G_{g}}$ are the shear modulus of the effective solid CNTs and graphite sheet, respectively, and ${\textstyle {\bar {J}}_{nt}}$ and ${\textstyle J_{ow}}$ are the polar moment of inertia of the effective solid CNTs and outer wall, respectively. They are

 ${\displaystyle {\bar {J}}_{nt}={\frac {\pi d_{nt}^{4}}{32}}\quad ,\quad J_{ow}={\frac {\pi \left(d_{nt}^{4}-\left(d_{nt}-2t\right)^{4}\right)}{32}}.}$
(4.24)

Replacing the expressions given by 4.24 in 4.23, the equivalent shear modulus can be written as

 ${\displaystyle {\bar {G}}_{nt}=\left[1-\left(1-{\frac {2t}{d_{nt}}}\right)^{4}\right]G_{g}.}$
(4.25)

Finally, it is necessary to obtain the new density of the effective solid CNTs, as the total weight of the MWCNTs can not change in the composite when they are considered a solid cylinder, then

 ${\displaystyle {\bar {\rho }}_{nt}={\frac {A_{nt}}{{\bar {A}}_{nt}}}\rho _{g}\quad \Rightarrow \quad {\bar {\rho }}_{nt}=\left[1-\left({\frac {d_{i}}{d_{nt}}}\right)^{2}\right]\rho _{g},}$
(4.26)

being ${\textstyle \rho _{g}}$ the density of the graphite sheet (${\textstyle \rho _{g}=2.25}$) and ${\textstyle d_{i}}$ the internal diameter of the MWCNTs.

The most common parameter used to define the amount of CNTs added to a composite is their weight fraction. However, the composite constitutive model developed requires knowing the volume fraction. The volume fraction of CNTs in the composite is the volume that occupies a solid cylinder with the same external diameter. This parameter can be calculated with the following expression [68] as

 ${\displaystyle k_{nt}={\frac {w_{nt}}{w_{nt}+{\frac {\bar {\rho _{nt}}}{\rho _{m}}}-{\frac {\bar {\rho _{nt}}}{\rho _{m}}}w_{nt}}},}$
(4.27)

where ${\textstyle w_{nt}}$ is the weight fraction and ${\textstyle \rho _{m}}$ is the density of the matrix.

### 4.2.4 Material non-linearity of the proposed model

In the proposed model, the composite performance is obtained from the mechanical response of its constituent materials, and each component is simulated with its own constitutive law. Therefore, if a constituent (i.e. the interface) is simulated with a non-linear law, the whole composite will become non-linear. As it has been already explained, with the present model it is possible to use any non-linear formulation to simulate the component behavior, such as plasticity, damage, viscosity, etc.

Besides the non-linear performance provided by each constituent, the load transfer capacity of the interface region is also affected if the interface is damaged. This effect must be included in the formulation.

According to Figure 3, the load is transferred from the interface to the CNTs reinforcement at their ends. Interface damage is expected to occur at the ends of the reinforcement, where there is larger stress concentrations. Assuming that the damaged region is unable to transfer loads and that the length required to transfer loads must remain constant, interface damage ends up affecting the parallel length of the nanotube, which can be calculated as

 ${\displaystyle l_{par}=l_{par}^{o}\left(1-d\right).}$
(4.28)

Here, ${\textstyle l_{par}^{o}}$ is the initial length of the nanotube working in parallel and ${\textstyle d}$ is the interface damage inner variable.

The dependence of the parallel length on the interface material damage provides a non-linear response of the composite, even when matrix and the carbon nanotube reinforcement are in their linear range.

## 4.3 Numerical Implementation

The proposed composite constitutive model has been implemented in PLCd [9], a finite element code that works with 3D solid geometries. The algorithm developed is described in Figure 4. PLCd has already implemented the constitutive laws that will be used to predict the performance of the composite components (elasto-plastic, elasto-damage and elastic). The formulation proposed has been written so that the constitutive laws of the constituents are seen as “black boxes”, following the recommendations of [33] and [4].

The FEM code enters into the new formulation with the prediction of the strain tensor of the composite material in the actual time step. Layers are assumed to have all the same strain; therefore the strain tensor of each layer is obtained rotating the composite strain to the direction in which the CNTs are oriented. In each layer, the strain of the matrix and the CNTs-interface are the same, as they work in parallel, (see 4.5). Knowing the strains for matrix material it is possible to obtain its stresses straightforward. On the other hand, to obtain the stresses for the CNT-interface material, it is necessary to separate it in two regions. In the flow chart shown in Figure 4, these two regions are represented as “Parallel Block” and “Serial Block”. This division is performed based on the value of ${\textstyle N^{par}}$ (defined in 4.12). This value depends on the damage evolution of the interface, as has been explained in section 4.2.4.

The Parallel Block corresponds to the central region, where the CNTs and the interface work in parallel behavior and, therefore, they have the same strains. In this region the stresses for each component are obtained from the strain tensor, using their constitutive equation. Finally, the stress tensor of the CNT-interface material in the “Parallel Block” at time ${\textstyle t+\Delta t}$ is

 ${\displaystyle \left[\sigma _{ntiz}^{par}\right]^{t+\Delta t}={\bar {k}}_{nt}\left[\sigma _{nt}^{par}\right]^{t+\Delta t}+{\bar {k}}_{iz}\left[\sigma _{iz}^{par}\right]^{t+\Delta t}.}$
(4.29)

On the other hand, at the ends of the CNTs, the interface-CNTs material has a serial behavior and it is necessary an initial prediction of the CNT or of the interface strains, in order to integrate the local stress in both components. If this initial prediction is made on the interface, its strains can be computed as

 ${\displaystyle \Delta \varepsilon _{ntiz}^{ser}=\left[\varepsilon _{ntiz}^{ser}\right]^{t+\Delta t}-\left[\varepsilon _{ntiz}^{ser}\right]^{t}}$
(4.30)
 ${\displaystyle \left[\Delta \varepsilon _{iz}^{ser}\right]_{o}=\left[C_{iz}^{ser}\right]^{-1}:C_{ntiz}^{ser}:\Delta \varepsilon _{ntiz}^{ser}}$
(4.31)
 ${\displaystyle \left[\varepsilon _{iz}^{ser}\right]_{o}^{t+\Delta t}=\left[\varepsilon _{iz}^{ser}\right]^{t}+\left[\Delta \varepsilon _{iz}^{ser}\right]_{o}.}$
(4.32)

And, the strain tensor of the interface in the iteration step ${\textstyle n}$ is used to calculate the strain tensor of the CNT as

 ${\displaystyle \varepsilon _{ntiz}^{ser}={\bar {k}}_{nt}\varepsilon _{nt}^{ser}+{\bar {k}}_{iz}\varepsilon _{iz}^{ser}}$
(4.33)
 ${\displaystyle \left[\varepsilon _{nt}^{ser}\right]_{n}={\frac {1}{{\bar {k}}_{nt}}}\left[\varepsilon _{ntiz}^{ser}\right]-{\frac {{\bar {k}}_{iz}}{{\bar {k}}_{nt}}}\left[\varepsilon _{iz}^{ser}\right]_{n}.}$
(4.34)

Once knowing the strain tensor of both component materials, the constitutive law of each one is used to calculate their stress tensor. Afterwards it is necessary to verify that the iso-stress condition is indeed fulfilled. Therefore

 ${\displaystyle \left[\Delta \sigma ^{ser}\right]_{n}=\left[\sigma _{iz}^{ser}\right]_{n}-\left[\sigma _{nt}^{ser}\right]_{n}\leq tolerance.}$
(4.35)

If the residual stress is greater than the tolerance, the prediction of the interface strain must be corrected. A Newton-Raphson scheme is adopted to do this correction. The method uses the Jacobian to update the unknown variable, in this case, the interface strain, then

 ${\displaystyle {\begin{array}{cl}J_{n}&=\left.{\frac {\partial \left[\Delta \sigma ^{ser}\right]_{n}}{\partial \varepsilon _{iz}^{ser}}}\right|_{\varepsilon _{iz}^{ser}=\left[\varepsilon _{iz}^{ser}\right]_{n}}={\frac {\partial \left[\sigma _{iz}^{ser}\right]_{n}}{\partial \varepsilon _{iz}^{ser}}}-{\frac {\partial \left[\sigma _{nt}^{ser}\right]_{n}}{\partial \varepsilon _{nt}^{ser}}}:{\frac {\partial \varepsilon _{nt}^{ser}}{\partial \varepsilon _{iz}^{ser}}}\\\\&=\left[C_{iz}^{ser}\right]_{n}-\left[C_{nt}^{ser}\right]_{n}\left(-{\frac {{\bar {k}}_{iz}}{{\bar {k}}_{nt}}}\right),\end{array}}}$
(4.36)

and, finally

 ${\displaystyle J_{n}=\left[C_{iz}^{ser}\right]_{n}+\left[C_{nt}^{ser}\right]_{n}\left({\frac {{\bar {k}}_{iz}}{{\bar {k}}_{nt}}}\right).}$
(4.37)

Therefore, the strain tensor of the interface for the next step n+1 is estimated as

 ${\displaystyle \left[\varepsilon _{iz}^{ser}\right]_{n+1}=\left[\varepsilon _{iz}^{ser}\right]_{n}-J_{n}^{-1}:\left[\Delta \sigma ^{ser}\right]_{n}.}$
(4.38)

This iterative process continues until the residual stress is smaller than the required tolerance.

The final stresses in the serial region “Serial Block” of the CNTs-interface are

 ${\displaystyle \left[\sigma _{ntiz}^{ser}\right]^{t+\Delta t}=\left[\sigma _{nt}^{ser}\right]^{t+\Delta t}=\left[\sigma _{iz}^{ser}\right]^{t+\Delta t}.}$
(4.39)

And at the end, the final stress tensor for a specific layer is obtained as

 ${\displaystyle {\begin{array}{cl}\left[\sigma \right]^{t+\Delta t}&=k_{m}\left[\sigma _{m}\right]^{t+\Delta t}+\\&\\&+\left(k_{nt}+k_{iz}\right)\left\{{\scriptstyle \left[N^{par}\right]^{t+\Delta t}\left[\sigma _{ntiz}^{par}\right]^{t+\Delta t}+\left[1-N^{par}\right]^{t+\Delta t}\left[\sigma _{ntiz}^{ser}\right]^{t+\Delta t}}\right\}.\end{array}}}$
(4.40)
 Figure 4: Flow chart of the proposed model in a FEM code.

# 5 Validation and numerical results

In this chapter, the validation of the proposed constitutive model using data from the literature is presented. Then, a numerical example is shown using the model calibrated. The basic formulation of the different constitutive models used for the simple materials can be seen in the Appendix A.

## 5.1 Validation of the elastic response

In the following section are compared the composite stiffness predicted by the composite constitutive model (see Section 4.2) with experimental data obtained from the literature. For this elastic properties validation the experimental data presented in the papers of Coleman et al. [65,70] is used. In these works several composites made of the same matrix with different MWCNTs are experimental tested.

#### Materials description

In the following, it will present the mechanical properties of the material components used and the composites data.

## Matrix component:

The matrix material is polyvinyl alcohol (PVA) and its Young's modulus is given by the authors as ${\textstyle E_{m}=1.9\pm {0.3}}$ [GPa] [65].

## Interface component:

The authors found that the Young's modulus of the crystalline polymer phase is of ${\textstyle E_{iz}=46}$ [GPa]. On the other hand, the parameter ${\textstyle {\frac {b}{r_{nt}}}}$ is estimated following the procedure described in Section 4.2.2.

## MWNTs component:

The nanotubes used in [65] are an arc grown MWCNT (Arc-MWCNT), two types of catalytic MWCNT from Nanocyl S.A. (CVD-1, CVD-2), a catalytic MWCNT produced in Orléans (France) (CVD-3), and a double walled nanotube (Dwnt). While in [70] the nanotube used is MWCNT from Nanocyl S.A. (MWCNT).

The maximum Young's modulus of the CNTs is ${\textstyle \sim {1}}$ [TPa][65], which corresponds to the stiffness of a perfect graphite sheet. The equivalent stiffness (see Section 4.2.3) of the nanotubes are calculated using this perfect stiffness value and considering a thickness of the outer layer of ${\textstyle t=0.34}$ [nm][34,68].

The most important collected data of the nanotubes used are presented in table 1:

 Type ${\displaystyle d_{nt}\left(\mathrm {nm} \right)}$ ${\displaystyle l_{nt}\left(\mathrm {\mu m} \right)}$ ${\displaystyle l_{nt}/d_{nt}}$ ${\displaystyle b/r_{nt}}$ ${\displaystyle {\bar {E}}_{nt}\left(\mathrm {GPa} \right)}$ ${\displaystyle N^{par}}$ Arc-MWCNT 24 1 42 0.81 56 0.97 CVD-3 16 3.8 238 1.47 83 0.99 CVD-2 14 2.1 150 2.27 95 0.99 CVD-1 15 1.8 120 2.83 89 0.98 Dwnt 2.5 2.2 880 4.87 470 0.99 MWCNT 15 1.72 115 3.30 89 0.98

## Composites:

A parameter missing in Table 1 is the direction distributions of the CNT. In general, obtaining this information from the composite is very complicated. To outstep this impediment it is possible to rewrite equation given by 4.11 for one layer as

 ${\displaystyle \mathbf {C} _{layer}=k_{m}\mathbf {C} _{m}+k_{ntiz}\mathbf {C} _{ntiz}^{eff},}$
(5.1)

where

 ${\displaystyle k_{ntiz}=k_{nt}+k_{iz}\qquad \mathbf {C} _{ntiz}^{eff}=N^{par}\mathbf {C} _{ntiz}^{par}+\left(1-N^{par}\right)\mathbf {C} _{ntiz}^{ser}.}$
(5.2)

Cox [71] and Krenchel [72] modified the rule of mixtures proposing the following equation to calculate the composite Young's modulus

 ${\displaystyle E=k_{m}E_{m}+k_{f}\eta _{o}E_{eff},}$
(5.3)

where ${\textstyle E_{m}}$ and ${\textstyle E_{eff}}$, are the Young's modulus of the matrix and effective reinforcement, respectively. The volume fraction for each component is ${\textstyle k}$ and ${\textstyle \eta _{o}}$ is a fiber orientation efficiency factor. For the present validation 5.3 will be modified, adapting it to the developed formulation. Therefore

 ${\displaystyle C_{composite}=k_{m}C_{m}+k_{ntiz}\eta _{o}C_{ntiz}^{eff}.}$
(5.4)

The value of the efficiency factor related to fiber orientation was taken from literature. In composites with a random distribution, ${\textstyle \eta _{o}=0.38}$.

#### Results

 Figure 5: Comparison of numerical and experimental results [65,70].

Figure 5 shows the values of ${\textstyle dC/dk_{nt}}$, this is: the slope of the curves of Young's modulus (${\textstyle C}$) divided by volume fractions of nanotubes (${\textstyle k_{nt}}$), for the different composites considered. In the figure the short lines represent the limits of the range experimental results presented in [65,70] and the red points correspond to the numerical result for each CNT type, obtained with the proposed composite model.

This figure shows that the formulation is capable of predicting the elastic stiffness of the composite, as most of the values obtained are comprehended between the limits defined by the experimental tests. There is only one case in which the value obtained exceeds the limits of the experimental test. This is because the effective Young's modulus of the Dwnt is highest since its diameter is really low.

## 5.2 Validation of the non-linear performance

The non-linear behavior of the composite constitutive model has been validated comparing the results provided by the model with the experimental data obtained from the paper of Meng et al. [73]. In that article the matrix used is Polyamide 6 (PA6) and all composites contained a 1 wt% of MWCNTs reinforcement.

The MWCNTs used in the experimental tests were purchased from Chengdu Organic Chemistry Co. Ltd. Two different composites where manufactured with these nanotubes. One of them contains the nanotubes “as is”, without any previous treatments. These nanotubes are called U-MWCNT. The other composite uses nanotubes that where treated with a mixture of concentrated sulfuric and nitric acids. These are called A-MWCNT.

#### Materials description

In the following, it will present the properties of the material components used and the information of the composites.

## Matrix component:

The matrix material is characterized with an isotropic, elasto-plastic model using a Von-Mises yield criterion. The mechanical parameters of the model were calibrated using the experimental data described in [73], obtaining a Young's modulus of ${\textstyle E_{m}=2.67}$ [GPa], a Poisson ratio of ${\textstyle \nu _{m}=0.4}$ and an elastic threshold of 35 [MPa]. The parameters used to simulate matrix material are validated comparing the stress-strain graph obtained with the numerical model with the experimental one. This comparison is shown in Figure 6.

 Figure 6: PA6 stress-strain relations for static tests [73].

## Interface component:

The interface zone is associated with the crystalline matrix around of MWCNTs. The properties of this material are better than those of the amorphous matrix. The volume fraction of the interface zone has been estimated with the data presented in the paper of Meng [73] and the equations developed in Section 4.2.2. On the other hand, the mechanical properties of the interface are used to calibrate the model. In current simulation, the interface has been defined with a isotropic, elasto-damage model with linear softening and Tresca yield surface. The mechanical parameters used are ${\textstyle E_{iz}=5}$ [GPa], ${\textstyle \nu _{iz}=0.4}$ and ${\textstyle G_{iz}=1.8}$ [GPa]. Damage in the interface starts for a stress threshold of 120 [MPa]. This value is in the range of theoretical and experimental tests value obtained in [46].

## MWCNTs component:

Numerical simulations of molecular structural mechanics of CNTs show that the Young's moduli are in the range of ${\textstyle 1.05\pm {0.05}}$ [TPa] and the shear moduli is about ${\textstyle 0.4\pm {0.05}}$ [TPa] [74]. It has been also shown that these values do not change significantly for CNTs with two, tree or four walls.

Regarding the transverse modulus of CNTs, it has been assessed from numerical and experimental results that there is an inverse relationship between axial and transverse modulus for carbon fibers [75]. Higher axial stiffness is associated to a longer and more aligned crystalline structure of the nanotube in this direction, which reduces properties in the transverse direction. Following this approach, in current simulation the transverse moduli of the MWCNTs are defined with the same values of the interface component.

Therefore, the equivalent properties of the MWCNTs were obtained using the equations described in Section 4.2.3. The diameter of MWCNT is ${\textstyle d_{nt}=50}$ [nm]. The measurement of several MWCNTs provided an estimation of the internal diameter of ${\textstyle d_{i}=8.2}$ [nm] [68]. The effective density of MWCNTs has a value of ${\textstyle {\bar {\rho }}_{nt}=2.2}$; and the volume fraction of MWCNTs in the composite is 0.51 % The MWCNTs have been simulated using an elastic orthotropic material with the following properties:

${\displaystyle E_{1nt}={\bar {E}}_{nt}=56}$[GPa] , ${\textstyle E_{2nt}=E_{3nt}=E_{iz}=5}$ [GPa]

${\displaystyle G_{12nt}=G_{13nt}={\bar {G}}_{nt}=41}$[GPa] , ${\textstyle G_{23nt}=G_{iz}=1.8}$ [GPa]

${\displaystyle \nu _{12nt}=\nu _{13nt}=\nu _{23nt}=\nu _{nt}=0.2}$

${\displaystyle \nu _{ij}={\frac {E_{i}}{E_{j}}}\nu _{ji}\quad \Rightarrow \quad \nu _{21nt}=\nu _{31nt}=0.018\quad \nu _{32nt}=0.2}$

## Composite:

The composites tested had a random distribution of the MWCNTs. This is simulated in the numerical model by dividing the composite in several layers, each one containing CNTs with a different orientation. Current simulation divides the composite in 10 layers and CNTs angles varying from ${\textstyle 0\,^{\circ }}$ to ${\textstyle 90\,^{\circ }}$. Each layer has a volume fraction of 10 % Table 2 shows the volume fractions of the three composite components in each layer. This table also shows some geometry information of the MWCNTs and the interface zone, as well as the initial value of ${\textstyle N^{par}}$.

 Composite ${\displaystyle k_{nt}}$[% ${\displaystyle k_{iz}}$[% ${\displaystyle k_{m}}$[% ${\displaystyle l_{nt}/d_{nt}}$ ${\displaystyle b/r_{nt}}$ ${\displaystyle N^{par}}$ PA6/A-MWCNT 0.5 4.1 95.4 250 2.00 0.98 PA6/U-MWCNT 0.5 5.3 94.2 250 2.35 0.98

#### Results

In Figure 7 are represented the numerical and experimental results obtained for the composite made with A-MWCNTs. This figure shows an initial reduction of the composite stiffness, result of matrix yielding. Afterwards damage begins in the interface zone and, consequently, the composite continues reducing its stiffness. Interface damage leads to a reduction of the parallel length (see 4.28). When the interface is completely damaged, the whole CNT-interface material has a serial performance. At this stage stresses in the interface are zero, and so must be the stresses in the carbon nanotubes. Therefore, the final stiffness of the composite corresponds to a material with a volumetric participation of 95.4 % of PA6 matrix, and the rest of the material correspond to voids.

 Figure 7: PA6/A-MWCNT stress-strain relations for static tests [73].

Figure 8 shows the results for the composite made with U-MWCNTs. This composite is the same than the previous one (made with A-MWCNTs), with the only difference that in this case the bond between U-MWCNTs and interface zone is weaker. To take into account this difference, the numerical model used for this composite is the same used for the previous one, varying the threshold at which damage starts in the interface. In current simulation this value is reduced to 70 [MPa].

This simulation provides a maximum stress in the composite lower than the value obtained for previous one, consequence of having a weaker interface. The simulation also shows some divergences between the numerical and the experimental values. Both graphs start to differ for a strain of 2.5 % and the maximum load reached by the numerical simulation is larger than in the experimental tests. However, it has to be noted that the experimental tests provide a maximum stress lower than having just plain matrix (see Figure 6). Therefore, the differences observed in Figure 8 may be justified.

 Figure 8: PA6/U-MWCNT stress-strain relations for static tests [73].

## 5.3 Simulation of a four points bending beam

### 5.3.1 Materials properties, geometry and FE model of the analyzed structure

This section presents the numerical simulation of a four points bending beam, which is shown in Figure 12. Although there is no experimental results available to compare the results obtained with the numerical simulations, the analysis performed is used to show the numerical performance of the proposed model and the effect of reinforcing the polymer with CNTs. The validity of the results is assumed, based on the comparisons made in Section 5.1 and 5.2. The composite used in the simulation is a reinforced matrix with MWCNTs and it has been proposed in the framework of M_RECT project. The matrix in the composite is a Polyether Ether Ketone (PEEK) thermoplastic polymer provided by Victrex company. While the MWCNTs are the NanocylNC7000 from Nanocyl S.A. Two composite with different weight fractions percentage (0.5 and 2.0 wt%) of MWCNTs reinforcement will be analyzed.

In the following, it will describe the properties and constitutive model of the material components and the composites information.

#### 5.3.1.1 Matrix component

An elasto-plastic constitutive model with hardening is applied to characterize the behavior of the PEEK component. The matrix material has a Young's modulus of ${\textstyle E_{m}=3.9}$ [GPA], a shear modulus of ${\textstyle G_{m}=1.9}$ [GPA], a Poisson ratio of ${\textstyle v_{m}=0.4}$. These elastic mechanical properties are obtained from M-RECT project and of the information provided by Victrex(http://www.victrex.com). The constitutive model is calibrated with an elastic threshold of ${\textstyle 32}$ [MPa] and an ultimate tensile strength of ${\textstyle 90}$ [MPa]. Figure 9 shows the comparison between the experimental data from the project and numerical results obtained with the constitutive model calibrated.

 (a) In a tensile test. (b) In a shear test. Figure 9: Comparison of the experimental data with the numerical results for PEEK.

#### 5.3.1.2 Interface component

The constitutive model used to simulate the behavior of the crystalline PEEK around of the MWCNTs is an elasto-damage model with exponential softening. The mechanic properties of this interface zone are obtained following the same procedure used by Coleman et al. [65] but using the information presented in the works of Diez-Pascual et al. [76,77], which use the PEEK material as matrix too. Then, the properties obtained are, a Young's modulus of ${\textstyle E_{iz}=5.07}$ [GPa], a shear modulus of ${\textstyle G_{iz}=2.47}$ [GPa] and a Poisson ratio of ${\textstyle v_{iz}=0.4}$. The value of the elastic threshold used in the model is of 28 [MPa]. This parameter is obtained when the composite constitutive model is calibrated to reproduce the experimental curve shown in Figure 10. This experimental data is for a PEEK reinforced with 3 wt% of MWCNTs obtained in the framework of the previous referred project.

 Figure 10: Experimental data and numerical response with the calibrated interface component model.

#### 5.3.1.3 MWCNTs component

The geometric characteristics of the MWCNTs are obtained from the paper of Jiang et al. [78], who obtains as average diameter and length of ${\textstyle 10.4}$ [nm] and ${\textstyle 0.7}$ [${\textstyle \mu }$m], respectively. For the simulation the MWCNTs are considered as an orthotropic elastic material. The equivalent properties are obtained using the equations described in Section 4.2.3, assuming ${\textstyle E_{g}=1.05\pm {0.05}}$ [TPa] and ${\textstyle G_{g}=0.4\pm {0.05}}$ [TPa][74] and a thickness of the outer layer of ${\textstyle t=0.34}$ [nm][34,68]. And taking the same consideration than before, the transverse properties are defined with the same values of the interface.

${\displaystyle E_{1nt}={\bar {E}}_{nt}=131}$ [GPa], ${\textstyle E_{2nt}=E_{3nt}=E_{iz}=5.07}$ [GPa],

${\displaystyle G_{12nt}=G_{13nt}={\bar {G}}_{nt}=104}$ [GPa], ${\textstyle G_{23nt}=G_{iz}=2.47}$ [GPa],

${\displaystyle \nu _{12nt}=\nu _{13nt}=\nu _{23nt}=\nu _{nt}=0.2}$

,

${\displaystyle \nu _{ij}={\frac {E_{i}}{E_{j}}}\nu _{ji}\quad \Rightarrow \quad \nu _{21nt}=\nu _{31nt}=0.008\quad \nu _{32nt}=0.2}$

.

#### 5.3.1.4 Composites

Table 3 shows the volume fractions of each component in the composites simulated.

 Composite ${\displaystyle k_{m}}$[% ${\displaystyle k_{nt}}$[% ${\displaystyle k_{iz}}$[% PEEK-0.5CNT ${\displaystyle 84.95}$ ${\displaystyle 0.35}$ ${\displaystyle 14.7}$ PEEK-2.0CNT ${\displaystyle 91.89}$ ${\displaystyle 1.41}$ ${\displaystyle 6.70}$

The orientation distribution of the MWCNTs has been defined assuming that the composite is formed by several layers, each one with a specific angle and volume fraction of MWCNTs. The volume fractions of the MWCNTs for the different layer in the composite are shown in the Figure 11. The value of the volume fractions in the figure are relative values respect to the total volume fraction of the MWCNTs in the composite.

 Figure 11: MWCNTs orientation distribution in the composite.

#### 5.3.1.5 Geometry of the simulated structure

The selected structure for the numerical simulations is a simple supported beam with two concentrated loads, which are applied at ${\textstyle 1/3}$ of both beam ends. Figure 12 shows the geometry and its dimensions, the boundary conditions and the load position on the analyzed structure.

 Figure 12: Geometry and extra information of the analyzed structure.

#### 5.3.1.6 FE model used

The symmetry of the geometry, of the applied load and of the boundary conditions of the structure allows to reduce the numerical model in the simulation. For this case, the reduced FE model is a quarter of the real geometry of the structure. Figure 13 shows the numerical model and the FE mesh used for the numerical analysis. The more relevant data about the FE mesh is shown in the Table 4.

 Item Nodes Elements Type Elem. Order Quantity/Type 1953 1200 Hexahedron Quadratic
 Figure 13: FE mesh used in the reduce model.

In order to obtain the real behavior of the structure with the reduced FE model it is necessary impose the restrictions on the numerical model given by the symmetry. There are two symmetry planes: The X-axis symmetry plane that has as normal axis the longitudinal axis (X-axis). In this symmetry plane, X direction displacements on the plane's nodes are restricted to zero. The other symmetry plane is the Y-axis symmetry plate, which has as normal axis the Y-axis in the model. For this symmetry plane, the null displacements restriction on the plane's nodes is Y direction. The longitudinal direction in the numerical model (X-axis) is taking as a reference direction for the definition of the layes' angle in the composite.

### 5.3.2 Linear analysis

The numerical results obtained in the simulation are presented in comparative form, taking as reference the result obtained for the non-reinforced matrix (plain PEEK material).

In all cases, the applied load for elastic analysis is a fixed displacement of a ${\textstyle -0.001}$ [mm] in Z direction at P position (see the Figure 12). The result considered for the comparison is the reaction force in Z direction on the support. As the imposed displacement is the same for all analysis, the reaction force increases when the material is the PEEK reinforced. Table 5 shows the results obtained for the different composites normalized by the non-reinforced PEEK results.

 PEEK PEEK-0.5%CNT PEEK-2.0%CNT 1 1.20 1.52

In the central section of the beam, between concentrated loads, there is a pure bending situation. While, at both ends of the beam there are a coupling bending and shear conditions. In Table 6, it is possible to observe the longitudinal (X direction) and shear (XZ direction) stresses obtained in the structure with the different composites.

 Longitudinal Stress Shear Stress PEEK PEEK-0.5%CNT PEEK-2.0%CNT

### 5.3.3 Non-linear analysis

In order to obtain the non-linear response of the structure the fixed displacement at P position is gradually increased in the simulation. Therefore, the reaction force in the support in Z direction increases too until the maximum value that the structure is able to take. The total force is four times the value obtained from the numerical model because the symmetry. The vertical (Z direction) displacement at the middle of the beam is taking as the reference loading increase.

The fixed displacement is applied to the numerical model in ${\textstyle 100}$ load steps of ${\textstyle 0.1}$ [mm]. Figure 14 shows the results obtained in the simulation for the different composites. When the vertical displacement is around ${\textstyle 1.5}$ [mm] the curves show the first loss of stiffness. This is because in the middle of the beam starts the plasticity in the PEEK. Then, when the vertical deformation is between ${\textstyle 3}$ [mm] to ${\textstyle 6}$ [mm] there is the second loss of stiffness. In this case, it is due to the damage in the interface zone. Subsequently, it is possible to observe that the model with non-reinforced PEEK has a higher stiffness than the ones with MWCNTs. This strange phenomenon is because at this point of the simulation the interface zone is completely damaged and, therefore, the contribution of the MWCNTs to the global stiffness is null. The stiffness obtained with the composites with MWCNTs are equivalent to the stiffness of a plain PEEK material but with some holes. This effect is clearly observed in Figure 15. This Figure shows the curves obtained for the simulation until a vertical displacement of ${\textstyle 50}$ [mm].

 Figure 14: No-linear structural response for PEEK-CNT.

Figure 15 shows a new loss of stiffness that takes place from ${\textstyle 30}$ [mm] to ${\textstyle 40}$ [mm] of vertical displacement. This laster loss of stiffness is because the plasticity model in the PEEK arrives to the ultimate tensile strength in the middle zone of the beam.

 Figure 15: No-linear structural response for PEEK-CNT up to ${\displaystyle 50}$ [mm] of vertical displacement.

Figure 16a shows the distribution of the longitudinal stress and Figure 16b of the shear stress for the composite reinforced with 0.5% of MWCNTs for a vertical displacement of 1 [mm]. The longitudinal plastic strain for the same composite and displacement state is shown in Figure 16c and the equivalent stress in Figure 16d.

 (a) Longitudinal stress distribution. (b) Shear stress distribution. (c) Longitudinal plastic strain distribution. (d) Equivalent stress distribution. Figure 16: No-linear results obtained at ${\displaystyle 1}$ [mm] of vertical displacement for PEEK-0.5%CNT.

### 5.3.4 Visco-elastic analysis

One of the main improvements shown by CNTs reinforced composites is their good damping response and energy dissipation, which makes them very useful for impact or vibration absorption purposes [69,79,80]. For this reason, in the following is analyzed the visco-elastic performance of the numerical model developed under such conditions. In order to obtain a viscous response of the composite it is necessary to use a visco-elastic model to characterize its components. The visco-elastic model used for the matrix and the interface zone is the generalized Maxwell model (see Figure 17) [3], which is already implemented in the PLCd code (see Section A.4 in Appendix A). To conduct the visco-elastic analysis, the MWCNTs are considered to have a linear elastic behavior.

 Figure 17: Scheme of the generalized visco-elastic Maxwell model [3].

 Mat./Prop. ${\displaystyle C_{_{\infty }}}$[MPa] ${\displaystyle C_{1}}$[MPa] ${\displaystyle \xi _{1}}$[MPa.seg] PEEK 3900 390 39 Interface 5070 1521 152

The response obtained, after calibrating the model (see Table 7), is shown in Figure 18. This figure shows the stress-strain curve of a point inside of the structure in a complete sinusoidal load-unload cycle. The numerical simulation has been conducted using the two composites previously described (0.5% and 2.0% wt) and the plain PEEK. Figure 18 shows that the areas enclosed by the curve, in the load-unload cycle, in the composites reinforced with MWCNTs are larger than the non-reinforced matrix (plain PEEK). In other words, the dissipation capacity of a composite with MWCNTs is better than the matrix alone. The composite reinforced with 0.5% wt of MWCNTs has higher dissipative capacity than the other composite (2.0% wt). This is because, in this composite, the volume fraction of the interface zone is higher than in the other, as it is shown in Table 3. This phenomenon occurs because when the volume fraction of the MWCNTs in the composite is low, their distribution in the composite improves, and then, increasing the interface volume.

## 5.4 Effect of the CNTs angle on the elastic properties

The present study about the influence of the CNTs angle in the resulting elastic properties using the developed composite constitutive model has been performed to asses the importance of having a good orientation of the CNTs in the composite.

The composite used for the analysis is a CNTs reinforced PEEK matrix and it has been also proposed in the framework of M_RECT project. The PEEK material used as matrix in the composite is the same than the one used in the above section (see Section 5.3.1). Therefore, the elastic properties previously obtained for the matrix and interface components have been used in this study. However, the MWCNTs used as reinforcements in the composite are the Baytubes C70P from Bayer, which have different geometric characteristics than the ones used in the above section. The current CNTs have an average diameter and length of ${\textstyle 13}$ [nm] and ${\textstyle 1}$ [${\textstyle \mu }$m], respectively. The CNTs are considered as an orthotropic elastic material and the equivalent properties are obtained using the equations described in Section 4.2.3, assuming ${\textstyle E_{g}=1.05\pm {0.05}}$ [TPa] and ${\textstyle G_{g}=0.4\pm {0.05}}$ [TPa][74] and a thickness of the outer layer of ${\textstyle t=0.34}$ [nm][34,68]. Considering that the transverse properties are defined with the same values of the interface, the CNTs properties are

${\displaystyle E_{1nt}={\bar {E}}_{nt}=105}$[GPa], ${\textstyle E_{2nt}=E_{3nt}=E_{iz}=5.07}$ [GPa],

${\displaystyle G_{12nt}=G_{13nt}={\bar {G}}_{nt}=85}$[GPa], ${\textstyle G_{23nt}=G_{iz}=2.47}$ [GPa],

${\displaystyle \nu _{12nt}=\nu _{13nt}=\nu _{23nt}=\nu _{nt}=0.2}$,

${\displaystyle \nu _{ij}={\frac {E_{i}}{E_{j}}}\nu _{ji}\quad \Rightarrow \quad \nu _{21nt}=\nu _{31nt}=0.008\quad \nu _{32nt}=0.2}$.

The composite simulated in the analysis has a 3% weight of CNTs, which represent a 1.94% of volume fraction. However, measurements made with X-rays show an apparent 5% weight. This difference is obtained because the MWCNTs have a higher apparent diameter than the original one. Therefore, the parameter ${\textstyle {\frac {b}{r_{nt}}}}$ is estimated assuming that this extra 2% weight in the measurement is the coating polymer around the nanotubes and then, it is the interface component considered. Taking this into account, the interface component has an approximate volume fraction of 1.31% and ${\textstyle {\frac {b}{r_{nt}}}=0.3}$.

 Figure 19: Change of elastic properties of CNTs reinforced PEEK with the CNTs angle.

Figure 19 shows the variation in the longitudinal elastic modulus, the transversal elastic modulus and the shear elastic modulus in the composite obtained when the CNTs angle is changed. The curves obtained, which are shown in the figure are symmetric with respect to the 45º vertical line. This behavior is different to the typical ones expected for reinforced matrices. In general, the maximum value of the longitudinal Young's modulus in the composite is obtained when the angle of the reinforcement is equal to 0º. For this reinforced PEEK with CNTs, the longitudinal Young's modulus increases from 0º until 40º and then, it decreases until its minimum value, which is obtained for an angle of 90º. A similar behavior presents the transversal Young's modulus but symmetrically respect to the 45º vertical line. The maximum value for the shear modulus is presented for an angle of 0º and the minimum for 45º.

# 6 Phenomenological homogenization. Concluding remarks

A new phenomenological composite constitutive model, based on the mixing theory, capable of predicting the mechanical performance of composites reinforced with carbon nanotubes has been presented. The formulation presented relates the reinforcement and the matrix in which they are embedded, using an interface material. This approach makes possible to consider non-linear phenomenons, such as debonding, by using non-linear constitutive laws to characterize the interface material. The formulation is written in a way in which all materials can be defined with their own constitutive law, improving the versatility of the model.

It has been shown that the elastic properties estimated with the model are in good agreement with experimental values obtained from literature. Only the numerical model of the composite made with the Dwnt reinforcement has given results in which the composite stiffness is overestimated. This is because the Dwnt has a very small diameter, which leads to a very high value of its equivalent Young's modulus.

The validation of the non-linear response provided by the proposed composite model has been performed using the experimental data of two different composites made with MWCNTs randomly distributed. The numerical curve obtained for the A-MWCNT is in good agreement with the experimental results. On the other hand, the numerical prediction obtained for the U-MWCNTs differs from the experimental results for strains larger than 2.5% However, it has to be said that the experimental results are lower than expected, as this composite is weaker than plain matrix.

Then, the formulation has used to predict and compare the mechanical properties of a straight beam subjected to four-point bending, with different material configurations. A non-linear analysis has also been made using the same structure and composites. The non-linear response of the beam obtained from the numerical simulation shows different points where there is a loss of structural stiffness. This structural behavior is obtained because the component materials in the composite reach their elastic thresholds and ultimate strength.

A visco-elastic material constitutive model is used for the polymeric matrix and the interface zone. The viscous response within the elastic range of the materials has been studied. The good capacity of energy dissipation in composites reinforced with MWCNTs has been proved with the simulations performed.

Finally, the composite constitutive model has been used to analyze the effect of nanotube orientation in the elastic properties of the composite, showing that the CNTs distribution within the reinforced matrix has a meaningful influence on the composite properties achieved.

All these tests have proved the validity of using a phenomenological model for the characterization of these materials. The developed composite constitutive model allows an accurate characterization of this kind of composites with an affordable computational cost, moreover taking into account the scale size of these reinforcements.

# 7 Introduction

In the last decades, several formulations have been developed to mathematically characterize and model composites as heterogeneous materials. Large number of composite models have been proposed to assess the global behavior of these materials by fulfilling the thermodynamic laws in linear and non-linear range. Constitutive equations have been developed for composites with different arrangements such as materials with long and short fibers, nanofibers, fibers laminated with one or more directions, and even random reinforcement distributions, etc. However, these formulations are limited because the constitutive relationships were made on a particular composite material and in general, can not be extrapolated to other composites.

The homogenization methods analyze the composite materials from an internal structure point of view. Over the years, many techniques have been developed, among them: the effective medium approach [81]; the self-consistent method [82]; the variational boundary method, which provides upper and lower limits of the total stiffness [83,84]; and the asymptotic homogenization method [85,86]. Due to the complex task that is required to represent the microscopic mechanical behavior of composites, homogenization approaches that use the RVE concept, together with stable computational methods, are very convenient in most cases. The behavior of the whole composite is obtained by a micromechanical study of the material components and their interaction within of the composite's microstructure through an RVE model.

Within the context of the homogenization methods, the known multiscale techniques use the RVE concept to address the characterization of composites and structures. As a result, these multiscale procedures do not obtain a closed form for the general constitutive equations. The stress-strain relationship of the composite is obtained by performing a detailed modeling of the its microstructure at the micro-scale in the RVE model. Among the main advantages of these techniques it is found that: they do not require any composite constitute assumption at the macro-scale; they can use any constitutive law in the simple materials, even non-linear response and time-dependency; and they can employ many computational techniques to find the response at the micro-scale, such as, the FEM, the Voronoi cell method, or numerical methods based on the fast Fourier transforms, etc.

The first-order homogenization [6] is one of the most extended and popular multiscale methods used nowadays. The approach uses the macro-scale deformation gradient tensor (or the strain tensor) to solve the micro-scale problem and then, by means of the microscopic results obtain the macro-scale stress tensor. The microscopic problem is solved through a Boundary Value Problem (BVP) on the RVE with particular boundary conditions, which are obtained using the macroscopic input. The resolution of the micro-scale BVP can be obtained through any mathematical or numerical approaches. After the solution of microstructural BVP, the microscopic displacement, deformation and stress fields are found and then, the macroscopic stress tensor is calculated as the volume average of the microscopic stress field. The BVP on the macro structure can be also approached through mathematical or numerical methods. When the solution of the coupled macro-micro BVPs is faced by FEM at both scales the formulation/implementation is called ${\textstyle {\hbox{FE}}^{2}}$ homogenization [7].

The advantages above described of the multiscale techniques become a challenge when a non-linear analysis is made on a realistic three-dimensional structure using a ${\textstyle {\hbox{FE}}^{2}}$ homogenization approach. In general, the required computational cost is extremely expensive for non-linear simulations because they require to solve, for each integration point at the macro-scale, one RVE at each time step of the analysis. Furthermore, from an energetic point of view the method must be consistent and therefore, the conservation of the dissipated energy through the scales should be guaranteed.

In the last decade, a second-order computational approach was proposed to be applied in critical regions of intense deformation, where the characteristic wave length of the macro-scale deformation field is of the order of the size of the micro-scale [8]. The homogenization technique has been developed as a natural extension of the first-order homogenization method. In the approach, the macroscopic gradient of the deformation gradient is also incorporated in the microscopic BVP and them, the stress tensor and a second-order stress tensor are retrieved. In the micro-scale level a first-order equilibrium problem is conserved but extra boundary conditions in the BVP must be verified. However, a full second gradient continuum theory appears in the macro-scale level, which requires solving a higher-order equilibrium problem. The solution of the reformulated macroscopic BVP is made through a complex finite element implementation, which has restricted its massive application.

The use of a second-order approach for the non-linear analysis of structures has the advantage of including higher-order effects on the micro-scale. Nevertheless, the complex computational implementation required to solve the macro-scale BVP restricts its application to realistic composite structures. From a mathematical and computational point of view the first-order approach is simpler than the second-order scheme. Then, it would be interesting to develop an enhanced-first-order approach which retains the easy computational implementation but with an enriched solution in the micro-scale.

## 7.1 Part's outline

In this second part of the present monograph a multiscale homogenization model for composite structures is addressed. The proposed formulation can simulate the behavior of three-dimensional structures in non-linear range.

Following this objective, in Chapter 8 a review of the state of the art of the more relevant homogenization theories is shown. The asymptotic homogenization theory introduces the concept of two or more scale lengths, therefore it lays the foundations of what is today known as multiscale homogenization using RVE concept. At the end of the chapter a review of the most important multiscale approaches found in the literature is made.

Chapter 9 presents the formulation of the two-scale homogenization developed in this study. In the beginning, the first-order homogenization approach is described and then, using concepts of the second-order approach, a enhanced-first-order homogenization is proposed. The BVPs in both scales and their computational implementation through of an efficient three-dimensional FEM is also shown.

In Chapter 10 the results obtained with the first-order homogenization proposed are validated with other microscopic formulations. From the comparison made, it can be observed that the multiscale formulation obtains a good characterization of the composite without any special consideration because the microstructure is modeled using an RVE. The computational run times and memory used are also compared showing the advantages and drawbacks of the proposed method.

Chapter 11 presents a novel strategy for the non-linear simulation of composite structures using multiscale homogenization approaches. This strategy is capable of reducing a 98% the computational time required in a non-linear analysis. The formulation is implemented in the FE code PLCd showing mesh independent at both scales as the develop approach conserves the energy through the scales.

Finally, in Chapter 12 the conclusions of the formulations presented in this part of the monograph are addressed in detail.

# 8 State of the art

The modeling of heterogeneous composite materials have a high degree of complexity because their constitutive behavior is strongly dependent on microstructural effects. The development of formulations based on multiple scales, capable of predicting the response of a heterogeneous material phenomenologically according to the information derived from a study at the micro-mechanical, is a natural choice to address the problem.

Numerous efforts have been made to mathematically model composite materials and structures with homogenization methods by using suitable multiscale techniques with relatively good approximation to the real global response of the composite. The most significant multiscale techniques, based in different theoretical principles, are:

## 8.1 The effective medium approximation

The method proposed originally by Eshelby [81] consists to find the current stresses in an elastic solid when a region of it, normally called inclusion (see Figure 20), suffers a change of shape and size which, if the surrounding material was absent, is represented by a uniform homogeneous deformation. By means of the Eshelby solution a number of very important boundary value problems can be solved, like: i) the solution for an ellipsoidal inclusion embedded within an elastically mismatched matrix, ii) the solution for an ellipsoidal cavity in an elastic solid and iii) the solutions for circular and elliptical cracks in an elastic solid. Moreover, several theories which estimate the elastic properties of composite materials use the Eshelby solution. Further developments by Hashin [87], obtain expressions for the elastic moduli and their values bounds of many-phase heterogeneous materials using an approximate method based on the variational theorems of the elasticity theory and on a concentric-spheres model. Besides, Mori and Tanaka [88] developed a method of calculating the average internal stress in the matrix of a material containing inclusions with transformation strain. The obtained actual stress in the matrix is the average stress plus the locally fluctuating stress, the average of which vanishes in the matrix. The developed method also shows that the average stress in the matrix is uniform throughout the material and independent of the position of the domain where the average treatment is carried out.

 Figure 20: Ellipsoidal region of the Eshelby inclusion [149].

## 8.2 The self-consistent method

It may be considered as an extension of the effective medium approximation. The method is proposed by Hill [82] and it uses similar concepts to the Hershey-Kröner theory of crystalline aggregates. The macroscopic elastic moduli of two-phase composites taking into account the inhomogeneity of stress and strain fields are estimated. It is required phases with the character of matrix and of effectively ellipsoidal inclusions, but they can have any concentrations in the composite. The model of three phases or generalized self-consistent model given by Christensen and Lo [89] is a more elaborate version. The method embeds the spherical (or cylindrical fibers) in a cap (spherical or cylindrical) which represents the elastic properties of the matrix. Then, this set is embedded in an infinite medium with the effective properties of the composite material to be determined. Finally, the effective elastic shear modulus of the material is obtained integrating the differential equations governing the behavior of the three-phase boundary conditions and of the applied loads.

A more recent work by Mercier and Molinari [90], which is based on the interaction law (postulated by Molinari [91] and validated by Mercier et al. [92] on the Eshelby problem) proposed two self-consistent schemes for perfectly disordered materials. The first one is valid for any non-linear behavior, and the second scheme is used to aggregates with phases having the same strain rate sensitivity. Both schemes predict accurately the overall response of the composite material and they are able to capture the strain and stress histories of the components too.

## 8.3 Bounding methods

The bounding methods provide the lower and upper limits to the total stiffness of the system or of the composites. In general, these methods obtain a simpler expressions for the effective elastic properties through the minimum potential and complementary energies. The most relevant of those are:

### 8.3.1 The classical bounds of Voigt and Reuss

The Voigt approach determines the elastic moduli by averaging stresses, expressed in terms of strains. The method assumes strain uniformity throughout the composite material. Therefore, the average strain of each phase is equal to the applied strain in the composite, which is similar to the rule of mixture assumption. Hence, considering that each component is linear-elastic, the following relationship can be obtained

 ${\displaystyle \mathbf {C} =\sum _{i=1}^{n}\mathbf {C} _{i}\cdot \phi _{i},}$
(8.1)

where ${\textstyle \mathbf {C} }$ is the effective elastic tensor of the composite, ${\textstyle \mathbf {C} _{i}}$ and ${\textstyle \phi _{i}}$ are the elastic constitutive tensor and the volume fraction of the ${\textstyle i-th}$ phase, respectively.

On the other hand, Reuss proposed to determine the elastic moduli by averaging strains, expressed in terms of stresses but assuming stress uniformity and then all components have the same stress. Therefore, the following estimation for the effective elastic tensor of the composite is obtained as

 ${\displaystyle \mathbf {C} =\left[\sum _{i=1}^{n}{\frac {\phi _{i}}{\mathbf {C} _{i}}}\right]^{-1}.}$
(8.2)

The Voigt and Reuss methods give an upper and lower bounds for the elastic moduli of a composite with an arbitrary random micro geometry. The Voigt approach gives the upper bound of the elastic moduli for the compound material, while the Reuss approach gives the lower bound. These bounding methods depend only of the phase volume fractions, and of course of the elastic constitutive tensors of the components, but do not require any further information or assumption in respect of the microstructure. The expressions obtained for the bounding effective moduli are valid even for anisotropic component materials

### 8.3.2 The variational bounding method

The authors Hashin and Shtrikman [83] use the variational formulation to obtained the upper and lower bounds for the effective elastic moduli of quasi-isotropic and quasi-homogeneous multiphase materials of arbitrary phase geometry. The obtained bounds are close enough to provide a good estimate to the effective moduli when the ratios between the different phase moduli are not too large. The method obtains analytical expressions for the elastic constants of a heterogeneous material with random isotropic distribution of phases. Later, Walpole [93,94] generalized the bounding method of Hashin and Shtrikman for materials with several phases, which may be arbitrarily anisotropic.

A variational method for bounding the effective properties of nonlinear composite materials with isotropic components with full variational principle status was proposed by Ponte Castañeda [84]. The author proposed two dual versions of the new variational principle and it is demonstrated their equivalence to each other and to the classical variational principles. The approaches are used to determine the bounds and to estimate the effective energy functions of nonlinear composites in the context of the deformation theory of plasticity. For completely anisotropic composites simpler forms of the classical bounds of Voigt and Reuss are recovered from the new variational principles. Also, for isotropic, particle-reinforced composites, as well as for transversely isotropic, fiber-reinforced composites simpler forms for nonlinear Hashin-Shtrikman bounds are obtained.

Lahellec and Suquet [95] proposed a new method for determining the overall behavior of composite materials which can be composed by nonlinear inelastic components. The evolution equations describing the constitutive behavior of the components can be reduced to the minimization of an incremental energy function by using an implicit time-discretization scheme. The alternative minimization problem is rigorously equivalent to a nonlinear thermoelastic problem with a transformation strain which is a nonuniform field. Comparisons with full-field simulations show that the present model is good as long as the variational procedure is accurate in the purely dissipative setting, when elastic deformations are neglected. If this is the case, the model accounts in a very satisfactory manner for the coupling between reversible and irreversible effects and is therefore an accurate model for treating nonlinear viscoelastic and elasto-viscoplastic materials. Lahellec and Suquet [96] in a second part, of the described work, proposed a proper modification of the second-order procedure of Ponte Castañeda and leads to replacing, at each time-step, the actual nonlinear viscoelastic composite by a linear viscoelastic one. The linearized problem is even further simplified by using an effective internal variable in each individual component. The resulting predictions are in good agreement with exact results and improve the predictions of the previous work. The analytical models presented in the previous sections are reasonably able to predict equivalent material properties for relatively simple geometries and low volume fraction of the second component inclusions. But they are, in general, incapable to obtain the evolution of stresses and strains in the microstructure. Moreover, the actual heterogeneous materials cannot be treated with these models because normally these composites have an arbitrary microstructural morphology.

## 8.4 The asymptotic homogenization theory

This theory has proven to be a powerful technique for the analysis of structural arrangements of two or more scale lengths. Through the use of asymptotic expansions of the displacement and stress fields, and appropriate variational principles, the homogenization methods can provide not only the effective (homogenized) material parameters, but also distributions of stresses and strains at the two levels. Bensoussan et al. [85] proposed an asymptotic expansion of the solution in terms of a parameter ${\textstyle \epsilon }$, which is the ratio of the period of the structure to a typical length in the domain. The link from the microscopic to the macroscopic description of the behavior of the system is given by solving the problem in two scales defined by the spatial variables ${\textstyle x}$ and ${\textstyle y}$, where ${\textstyle x}$ is a macroscopic quantity and ${\textstyle y=x/\epsilon }$ is a microscopic quantity; ${\textstyle y}$ is associated with the small length scale of the inclusions or heterogeneities. The asymptotic problem is formulated in mathematical terms as a family of partial differential operators, depending on the small parameter ${\textstyle \epsilon }$. The operators may be time independent or time dependent, steady or of evolution type, linear or nonlinear, etc. The coefficients of the operators are periodic functions in all or in some variables with periods proportional to ${\textstyle \epsilon }$. Since ${\textstyle \epsilon }$ is assumed to be small, it is has a family of operators with rapidly oscillation coefficients. The standard classical formulation of this theory is found in the work of Sanchez-Palencia [86,97]. The two-scale process introduced in the partial differential equations of the problem produces equations in ${\textstyle x}$, in ${\textstyle y}$ and in both variables. Normally, the equations in ${\textstyle y}$ are solvable if the microscopic structure is periodic, and this leads to deduction of the macroscopic equations for the global behavior in ${\textstyle x}$. It is emphasized that the homogenized coefficients depend on the local structure of the medium. Fish et al. [98] presented a generalization of the mathematical homogenization method based on double-scale asymptotic expansion to account for damage effects in heterogeneous media. A closed-form expression relating local fields to the overall strain and damage is derived. Non-local damage theory is developed by introducing the concept of non-local phase fields (stress, strain, free energy density, damage release rate, etc.) in a manner analogous to that currently practiced in concrete. Numerical results of the model were found to be in good agreement with experimental data.

## 8.5 Homogenization using the RVE concept

The multiscale homogenization based on the use of a unit cell or RVE has emerged as one of the most promising methods to compute the response of composite structures. The unit cell is defined as a microscopic subregion that is representative of the entire microstructure in an average sense. The RVE is employed to obtain the effective properties for the homogenized material because it is assumed that it must contain a sufficient number of heterogeneities.

A thorough examination of the several representative volume element definitions found in the literature can be consulted in the review made by Gitman [99]. Moreover, Ostoja-Starzewki [100] points out that the RVE is very clearly defined in two situations only: i) unit cell in a periodic microstructure, as shown Figure 21a and ii) volume containing a very large (mathematically infinite) set of sub-scale elements, possessing statistically homogeneous and ergodic properties as shown Figure 21b.

 Figure 21: Representative volume element models: (a) unit cell approach; (b) statistical and ergodic approach [149].

On the other hand, Suquet [6] gives the basic principles of homogenization to obtain the constitutive equations for homogenized properties of a heterogeneous material. The process is resumed in the following three-step scheme and in Figure 22.

• Definition of a representative volume element. The size of this should be large enough to contain a sufficient number of micro heterogeneities of which the constitutive behavior on these individual constituents is assumed to be known.
• Microscopic boundary conditions are obtained based on the macroscopic input variables (e.g. strain tensor), taking into account the geometry, constitutive laws, etc.
• Macroscopic output variables are obtained based on the solution of the microscopic behavior of the RVE. Macroscopic properties of the equivalent homogeneous medium are evaluated.
 Figure 22: Representative homogenization scheme [149].

Terada et al. [101] conducted a study for multiscale homogenization of the convergence of the overall material properties when the unit cell size is increased. For general heterogeneous media that reveal irregular material distribution in a micro scale, the study showed that the periodic boundary condition provides the most reasonable estimates among the class of possible boundary conditions for statistical homogeneous media. The work demonstrated that it is not strictly required the periodicity of RVE geometry in evaluating the effective properties. Moreover, at non-linear range the mechanical behaviors are more sensitive to the size of RVE than those of linear case. The authors concluded that with more accurate geometric models of the RVE, and with larger RVE regions analyzed, the analysis provides a better understanding of the actual phenomena in the microscopic region, as well as a better characterization of the overall mechanical material performance.

### 8.5.1 Multiscale computational homogenization

Renard and Marmonier [102] were the first to use a finite element discretization to model heterogeneous materials with a multiscale approach. The method consists in solving two finite element problems, one for each scale. In the micro scale the geometry of the RVE is meshed and homogenization rules are used to link this with the macro scale problem.

Guedes and Kikuchi [103] studied the mechanical behavior of linear elastic 2D and 3D composite materials through the homogenization method. This study was one of the first that used homogenization in three dimensions. The method was implemented via a finite element technique to analyze a composite with a periodic microstructure. The implementation enabled the calculation of homogenized material mechanical properties, these characterized the overall behavior of the composite, as well as the analysis of the local performance of the material, as it enables the computation of local stress and strain distribution within the microstructure of the composite. An adaptive finite element method was introduced in order to improve the accuracy of the numerical results. An error measure is suggested for the homogenized material constants, based on the a priori error estimations and the numerical implementation.

#### 8.5.1.1 Non-linearity and FE2 homogenization

Swan [104] presented a computational stress and strain controlled homogenization methods for inelastic periodic composites within the framework of displacement finite element. The implementation of the stress controlled method employs a penalty formulation to insure that the displacement solution satisfies the linear-periodic decomposition. The methods were assessed on a complete set of homogenization computations for an elasto-plastic composite, the strain-controlled method was easier to implement and was demonstrated that it has more computationally efficient than the others.

Later, Smit et al. [7] presented a homogenization method for large deformations and viscoelastic material behavior on microscopic and macroscopic level. The homogenization method was implemented in a multi-level two dimensional finite element program with meshes on macroscopic level (mesh of entire structure) and microscopic level (meshes of RVEs). The local macroscopic stress is obtained by applying the local macroscopic deformation on a unique RVE through imposing appropriate boundary conditions and averaging the resulting stress field. To each macroscopic integration point a unique discretized RVE is assigned that provides the local macroscopic stress tensor and the tangential stiffness matrix. A separate iterative finite element procedures on the RVE is used in each iteration cycle of each macroscopic increment. The following Figure 23 shows the proposed scheme by the authors. The computational cost was the main disadvantage of this implementation. The increase in computational time with respect to a single-scale analysis on the same macroscopic mesh was of 160 times for the structure analyzed.

 Figure 23: Schematic representation of a multiscale finite element program [7].

Michel et al. [105] presented a review of several problems which are specific of composites with periodic microstructure composed of linear or non-linear constituents. The study obtains the estimation of phenomenological macroscopic constitutive models through the analysis of a microscopic representative volume element proposing the concept of “macroscopic degrees of freedom”. A general framework permitting either a strain or stress control was proposed, and the implementation of different types of boundary conditions was presented.

Haj-Ali and Muliana [106] proposed a three dimensional micromechanical modeling approach for the non-linear viscoelastic behavior of pultruded composites. A sub-laminate model is used to provide for a nonlinear equivalent continuum of a layered medium. The system is idealized using a weighted-average response of two simplified micromodels with fiber and matrix. The proposed micromechanical framework was able to generate the effective anisotropic non-linear viscoelastic response as a direct outcome of the micromechanical homogenization subcell. The same authors [107] extended the method to an integrated micromechanical and structural framework for the non-linear viscoelastic analysis of laminated composite materials and structures. The micromechanical model is numerically implemented within a shell-based non-linear finite element by imposing a plane stress constraint on its 3D formulation. The micromechanical model provides the effective non-linear constitutive behavior for each Gauss point. The formulation was validated with several experimental creep tests available in the literature. Finally, the work presents examples for non-linear viscoelastic structures, like a laminated panel and a composite ring. In a latter work, Haj-Ali et al. [108] extend the above described method for the analysis of thick-section fiber reinforced plastic (FRP) composite materials and structures. The proposed modeling framework is applied to a pultruded composite system. Non-linear 3D micromechanical models representing the different composite layers are used to generate through-thickness composite's effective responses. The nested non-linear micromechanical models are implemented at each integration point in the finite element structural analysis. The results obtained demonstrated good prediction capabilities for effective properties and for multi-axial non-linear behavior of pultruded composites.

Matsui et al. [109] made a feasibility study and introduced a parallel algorithm to achieve the computational efficiency. The focus was to the inhomogeneous deformation of the overall structure, which may imply the loss of periodicity assumed in the initial state. The study presented an efficient algorithm for the deconcentration of computational loads by using a PC-cluster system. A simple numerical example for three dimensional heterogeneous structure was made. The authors concluded that the two-scale analysis is still expensive, even if a PC-cluster system for parallel computations is available at non-linear range.

Ladeveze et al. [110,111] proposed a mixed, multilevel domain decomposition method or more accurately, as a “structure decomposition” method. The multiscale computational strategy consists of describing the structure as an assembly of simple components: substructures and interfaces. Each of these entities has its own variables and equations. The distinction between the micro and macro levels is made only at the interfaces, where forces and displacements are split into macro contributions and micro complements. A substructure is subjected to the action of its environment (the neighboring interfaces) defined by the force and a velocity distribution on its boundary. An interface between two substructures transfers both the velocity and the force distributions.

Oller et al. [112] extended the work presented by Zalamea [113], with a two-scale numerical homogenization method that assumes the periodicity of the internal structure of the material to address the problem of the steep gradient in the macroscopic field. To prevent the steep gradient in the macroscopic field variables, which is produced by local boundary effects, the authors proposed a local refinement of the finite element mesh. The objective was to maintain the periodic condition on the boundaries of the cells near to the perturbation.

De Souza Neto and Feijoo [114] discussed some equivalence relationships for large strain multiscale solid constitutive models based on the volume averaging of the microscopic stress and deformation gradient fields over a representative volume element. The work established that the volume averaging of the first Piola-Kirchhoff stress over the reference configuration of an RVE is mechanically equivalent to the spatial averaging of Cauchy stress over the deformed RVE configuration. Whenever such conditions are met, multiscale constitutive models resulting from either the reference or spatial stress averaging are identical.

#### 8.5.1.2 High-order multiscale homogenization

The multiscale methods mentioned in the above sections are commonly denominated as first-order homogenization approaches because in the mathematical formulation considers only the first gradient of the macroscopic displacement field.

Kouznetsova and Geers proposed what is called second-order homogenization [115,8], which is an extension of the first-order theory. In this case, the macroscopic deformation gradient tensor and its Lagrangian gradient is used to solve the boundary value problem at the microstructural scale. The second-order approach allows solve problems in the presence of localization phenomena without loss of precision in the solution because the Lagrangian tensor is taken into account. The main drawbacks of this method are its computational cost and complex implementation.

Geers, Coenen and Kouznetsova [116,117] proposed a computational homogenization technique for thin-structured sheets based on the homogenization for first and second-order continua. The three dimensional heterogeneous sheet is represented by a homogenized shell continuum for which the constitutive response is obtained from the nested analysis of a microstructural representative volume element, incorporating the full thickness of the sheet and an in-plane representative cell of the macroscopic structure. At an in-plane integration point of the macroscopic shell, the generalized strains (the membrane deformation and the curvature) are used to formulate the boundary conditions for the micro RVE problem. All microstructural constituents are modeled at the RVE as an ordinary 3D continuum. Upon proper averaging of the RVE response, the macroscopic generalized stress and the moment resultants are obtained.

# 9 Multiscale homogenization formulations

## 9.1 Introduction

In the context of solid mechanics and multiscale computational homogenization, one of the most extended and popular method is called first-order computational homogenization [105,101,118]. In this approach, the macro-scale strain tensor (or deformation gradient tensor) is used as input to solve the micro-scale Boundary Value Problem (BVP). The material stress-strain relationship is obtained from the solution of problem at the micro-scale i.e. the RVE which contains the detailed modeling of the internal heterogeneous structure of the composite. Therefore, it does not require any composite constitutive assumption or compatibility equation to address the composite response [119]. And, there are not restriction on the constitutive model used in the component materials, even non-linear materials and time-dependency models can be taken into account [7,112,120].

In the last decade, a second-order computational homogenization was proposed as a natural extension of the first-order approach [115,8,121]. It was developed to be applied in critical regions of large gradient deformation, where the characteristic wave length of the macro-scale deformation field is of the order of the size of the micro-scale. In this method the macroscopic gradient of the deformation gradient is also incorporated as input in the micro-scale BVP. The first-order equilibrium problem is conserved at the micro-scale level, while a higher-order equilibrium problem appears at the structural scale. The finite element framework necessary for the numerical solution of the macro-scale problem leads to many complications [122], which has restricted its extensive applicability.

The second-order homogenization approach is able to capture the second-order effects in the microscopic scale due to macroscopic high-order phenomena such as bending or strain localization, this is its major improvement over the first-order approach. However, the first-order homogenization conserves first-order equilibrium equations at both scales, which represents an advantage from a computational point of view. With the purpose to take the best of these two methods, the presented work proposes an enhanced-first-order homogenization approach [123]. The developed procedure takes into account macroscopic second-order effects by using the macro-scale second-order deformation measure in the micro-scale BVP. Besides, the proposed formulation conserves the classical first-order equilibrium problem in the structural scale.

## 9.2 General considerations

The no lineal transformation between the reference configuration of the body ${\textstyle \Omega }$ and the current configuration of the same body ${\textstyle \Omega ^{c}}$ is defined as: ${\textstyle {\boldsymbol {\phi }}\colon \Omega \rightarrow \Omega ^{c}\mid \mathbf {x} =\mathbf {\boldsymbol {\phi }} \left(\mathbf {X} \right)}$, where ${\textstyle \mathbf {x} \in \Omega ^{c}}$ and ${\textstyle \mathbf {X} \in \Omega }$ are respectively the current and the reference positions of the material point. Therefore, the linear mapping for an infinitesimal material line element is

 ${\displaystyle d\mathbf {x} =\mathbf {F} \cdot d\mathbf {X} ,}$
(9.1)

where the deformation gradient tensor is defined by

 ${\displaystyle \mathbf {F} ={\frac {\partial {\boldsymbol {\phi }}}{\partial \mathbf {X} }}=\nabla \mathbf {x} .}$
(9.2)

Here, the gradient operator ${\textstyle \nabla }$ is taken respect to the reference configuration ${\textstyle \mathbf {X} }$.

Nevertheless, if now a finite material line within a finite volume is considered, the expression given by 9.1 does not apply any more. However, a Taylor series expansion (centered at ${\textstyle \mathbf {X} _{o}}$) can be used to obtain an expression for the finite material line ${\textstyle \triangle \mathbf {x} }$ in the current configuration as

 ${\displaystyle \triangle \mathbf {x} =\mathbf {F} \left(\mathbf {X} _{o}\right)\cdot \triangle \mathbf {X} +{\frac {1}{2}}\mathbf {G} \left(\mathbf {X} _{o}\right):\triangle \mathbf {X} \otimes \triangle \mathbf {X} +{\mathcal {O}}\left(\triangle \mathbf {X} _{o}^{3}\right),}$
(9.3)

 ${\displaystyle \mathbf {G} ={\frac {\partial }{\partial \mathbf {X} }}\left({\frac {\partial {\boldsymbol {\phi }}}{\partial \mathbf {X} }}\right)=\nabla \mathbf {F} .}$
(9.4)

It can be shown that the third-order tensor ${\textstyle \mathbf {G} }$ has the symmetry property ${\textstyle G_{ijk}=G_{ikj}}$.

## 9.3 First-order homogenization approach

Let us consider a solid domain (or body ${\textstyle \Omega }$) with a periodic or quasi-periodic microstructure that can be represented by a Represent Volume Element. In this body, it is possible to establish two scale levels, a macro scale (or structural scale) for the macrostructure, and the other one micro scale (or sub scale) for the microstructure. The microstructural scale is defined using a RVE which characterizes the microstructure of the material. Let us also consider an infinitesimal material point ${\textstyle \mathbf {X} _{o}}$ in the reference configuration of the structure, and the RVE around this considered point as Figure 24 is showing.

 Figure 24: Macrostructure and microstructure around of the point ${\displaystyle \mathbf {X} _{o}}$.

The called principle of separation of scales [124] establishes that: the microstructural length scale ${\textstyle l_{\mu }}$ is assumed to be much smaller than the macrostructural characteristic length ${\textstyle l}$, which is the length over the macroscopic space. In other words, the principle says that the existing periodical microscopic dimension around of the macrostructural point (${\textstyle \mathbf {X} _{o}}$) must be smaller than the characteristic macrostructural dimension. If this principle is satisfied, the current configuration or deformed position of a material point in the RVE ${\textstyle \mathbf {x} _{\mu }\in \Omega _{\mu }^{c}}$ can be approximated as

 ${\displaystyle \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {x} _{\mu }^{o}+\mathbf {F} \left(\mathbf {X} _{o}\right)\cdot \triangle \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right),}$
(9.5)

where ${\textstyle \triangle \mathbf {X} _{\mu }=\mathbf {X} _{\mu }-\mathbf {X} _{\mu }^{o}}$, and ${\textstyle \mathbf {X} _{\mu }\in \Omega _{\mu }}$ is the reference configuration or non-deformed position of the material point in the RVE and ${\textstyle \mathbf {X} _{\mu }^{o}}$ and ${\textstyle \mathbf {x} _{\mu }^{o}}$ are the origin of the reference and the current coordinate system on the RVE, respectively (see Figure 25). The extra term ${\textstyle \mathbf {w} }$ is a microstructural displacement fluctuation field.

To simplify the symbolic manipulation of the formulation is convenient to set the coordinate system's origin as

 ${\displaystyle \mathbf {X} _{\mu }^{o}=0\qquad and\qquad \mathbf {x} _{\mu }^{o}=0.}$
(9.6)

Later, it will be proved that with these values, the rigid body motion of the RVE is avoided. Considering these restrictions, the expression given by 9.5 can be rewritten as

 ${\displaystyle \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {F} \left(\mathbf {X} _{o}\right)\cdot \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right).}$
(9.7)
 Figure 25: Reference and current configuration of the RVE.

### 9.3.1 Displacement field on the RVE

The displacement field ${\textstyle \mathbf {u} _{\mu }}$ at the RVE is defined by

 ${\displaystyle \mathbf {u} _{\mu }=\mathbf {x} _{\mu }-\mathbf {X} _{\mu },}$
(9.8)

and taking into account 9.7 in the previous equation,

 ${\displaystyle \mathbf {u} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \left[\mathbf {F} \left(\mathbf {X} _{o}\right)-\mathbf {I} \right]\cdot \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right),}$
(9.9)

where ${\textstyle \mathbf {I} }$ is the second-order unit tensor.

### 9.3.2 Kinematically admissible displacement fields and boundary conditions in the RVE

The displacement fields in the RVE that are kinematically admissible are obtained as a result of the coupling between the macrostructure and the microstructure. This linkage is based on the average theorems and they have been initially proposed for infinitesimal deformations by Hill [5]. Later, Hill [125] and Nemat-Nasser [126] extended these to finite deformations.

The first of the averaging relations postulated that the volume average of the microstructural deformation gradient tensor ${\textstyle \mathbf {F} _{\mu }}$ over the RVE must be equal to the macroscopic deformation gradient tensor ${\textstyle \mathbf {F} }$. In the considered point ${\textstyle \mathbf {X} _{o}}$ this is

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV,}$
(9.10)

where ${\textstyle V_{\mu }}$ is the volume of the RVE in the reference configuration.

Considering 9.7 it is possible to obtain the deformation gradient tensor in the microstructural scale as

 ${\displaystyle \mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)=\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {F} \left(\mathbf {X} _{o}\right)+\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right),}$
(9.11)

and using this relation, the right hand size of 9.10 is

 ${\displaystyle {\begin{array}{rcl}\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV&=&\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV,\\[6mm]&=&\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)+{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV.\end{array}}}$
(9.12)

Equation 9.12 can be rewritten as

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV,}$
(9.13)

or

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV.}$
(9.14)

Finally, applying the divergence theorem, in the right hand size of 9.14, this can be also rewritten in term of surface integral as

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA-{\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA,}$
(9.15)

where ${\textstyle \partial \Omega _{\mu }}$ is the RVE boundary domain in the reference configuration, and ${\textstyle \mathbf {N} }$ denotes the outward unit normal on ${\textstyle \partial \Omega _{\mu }}$.

Clearly, to satisfy the first average theorem, the integrals that depend of the displacement fluctuation in both 9.14 and 9.15 must vanish. Therefore,

 ${\displaystyle \int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV=\mathbf {0} }$
(9.16)

and

 ${\displaystyle \int _{\partial \Omega _{\mu }}\mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA=\mathbf {0} .}$
(9.17)
 Figure 26: Normal vectors to the surfaces in the reference configuration of a Cubic RVE.

Noting Figure 25 and considering that the reference geometry configuration of the RVE is originally a cube, as the figure is showing, the integral restriction on the RVE boundary can be splitted in the different surfaces of the ${\textstyle \partial \Omega _{\mu }}$ domain. Besides, taking the reference coordinate system that is shown in Figure 26, the outward unit normal of the cubic faces satisfy: ${\textstyle \mathbf {N} _{X}^{-}=-\mathbf {N} _{X}^{+}}$, ${\textstyle \mathbf {N} _{Y}^{-}=-\mathbf {N} _{Y}^{+}}$ and ${\textstyle \mathbf {N} _{Z}^{-}=-\mathbf {N} _{Z}^{+}}$. Here, the subscript makes reference to the axis which is perpendicular to the considered face and the superscript defines the position of the face on the axis. Therefore, considering this geometry, the expression given by 9.17 may be rewritten as

 ${\displaystyle {\begin{array}{l}(\int _{\mathbf {N} _{X}^{+}}\mathbf {w} ~dA_{yz}-\int _{\mathbf {N} _{X}^{-}}\mathbf {w} ~dA_{yz})\otimes \mathbf {N} _{X}^{+}\\[5mm]+(\int _{\mathbf {N} _{Y}^{+}}\mathbf {w} ~dA_{xz}-\int _{\mathbf {N} _{Y}^{-}}\mathbf {w} ~dA_{xz})\otimes \mathbf {N} _{Y}^{+}\\[5mm]+\int _{\mathbf {N} _{Z}^{+}}\mathbf {w} ~dA_{xy}-\int _{\mathbf {N} _{Z}^{-}}\mathbf {w} ~dA_{xy})\otimes \mathbf {N} _{Z}^{+}=\mathbf {0} .\end{array}}}$
(9.18)

Equation 9.18 shows that the boundary restriction on the displacement fluctuation field can be splitted on the different surface pairs (X, Y and Z) of the RVE boundary.

Previous equations from 9.16 to 9.18 can be used to obtain the different displacement fluctuation fields kinematically admissible in the microstructural level. Several models have been defined that assume different fluctuation fields, they are mentioned in the following.

(i) Taylor model (or zero fluctuations): The expression given by 9.16 is verified when

 ${\displaystyle \mathbf {w} ,\;{\hbox{sufficiently regular}}\vert \mathbf {w} \left(\mathbf {X} _{\mu }\right)=\mathbf {0} ,\;\forall \;\mathbf {X} _{\mu }\in \Omega _{\mu }.}$
(9.19)

This model gives homogeneous deformation in the microstructural scale level (see 9.24).

(ii) Linear boundary displacements (or zero boundary fluctuations): The expression given by 9.17 is verified when

 ${\displaystyle \mathbf {w} ,\;{\hbox{sufficiently regular}}\vert \mathbf {w} \left(\mathbf {X} _{\mu }\right)=\mathbf {0} ,\;\forall \;\mathbf {X} _{\mu }\in \partial \Omega _{\mu }.}$
(9.20)

The deformation of the RVE boundary domain for this class are fully prescribed.

(iii) Periodic boundary fluctuations:

The key kinematical constraint for this class is that the displacement fluctuation must be periodic on the different faces of the RVE. That is, for each pair ${\textstyle \left\{\mathbf {X} _{\mu }^{+},\mathbf {X} _{\mu }^{-}\right\}}$ of boundary points the expression given by 9.18 is verified when

 ${\displaystyle \mathbf {w} ,\;{\hbox{sufficiently regular}}\vert \mathbf {w} \left(\mathbf {X} _{\mu }^{+}\right)=\mathbf {w} \left(\mathbf {X} _{\mu }^{-}\right),\;\forall \;pairs\;\left\{\mathbf {X} _{\mu }^{+},\mathbf {X} _{\mu }^{-}\right\}\in \partial \Omega _{\mu }.}$
(9.21)

(iv) Minimal constraint (or uniform boundary traction):

In this constraint the nontrivial solution of 9.17 is obtained.

### 9.3.3 Microscopic and macroscopic strain tensor

Considering a infinitesimal deformation framework the strain tensor in the microstructural level can be obtained as

 ${\displaystyle {\begin{array}{rcl}\mathbf {E} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)&=&{\frac {1}{2}}\left(\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)+\mathbf {F} _{\mu }^{T}\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\right)-\mathbf {I} \\[3mm]&=&{\frac {1}{2}}\left(\mathbf {F} \left(\mathbf {X} _{o}\right)+\mathbf {F} ^{T}\left(\mathbf {X} _{o}\right)\right)-\mathbf {I} \\[3mm]&+&{\frac {1}{2}}\left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)+\left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\right)^{T}\right),\end{array}}}$
(9.22)

and, if 9.10 is satisfied it can be proved that taking the volume average of the microscopic strain tensor over the RVE domain the following relationship is obtained,

 ${\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {E} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV={\frac {1}{2}}\left(\mathbf {F} \left(\mathbf {X} _{o}\right)+\mathbf {F} ^{T}\left(\mathbf {X} _{o}\right)\right)-\mathbf {I} =\mathbf {E} \left(\mathbf {X} _{o}\right).}$
(9.23)

Here, ${\textstyle \mathbf {E} \left(\mathbf {X} _{o}\right)}$ is the macroscopic strain tensor. Therefore, it is possible to rewrite 9.22 as

 ${\displaystyle \mathbf {E} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)=\mathbf {E} \left(\mathbf {X} _{o}\right)+\mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right),}$
(9.24)

where ${\textstyle \mathbf {E} _{\mu }^{w}={\frac {1}{2}}\left(\nabla \mathbf {w} +\left(\nabla \mathbf {w} \right)^{T}\right)=\nabla ^{s}\mathbf {w} }$ is the contribution of the displacement fluctuation field to the microscopic strain tensor and ${\textstyle \nabla ^{s}}$ is the symmetric gradient operator. Because 9.10 is verified the volume average of ${\textstyle \mathbf {E} _{\mu }^{w}}$ over the RVE domain is equal zero.

### 9.3.4 Hill-Mandel principle and RVE equilibrium

The Hill-Mandel energy condition [82,127], also referred to as the macro-homogeneity condition, states that the virtual work of the point ${\textstyle \mathbf {X} _{o}}$ considered must be equal to the volume average of the virtual work in the RVE to any kinematically admissible displacement field, this principle can be formulated as

 ${\displaystyle \mathbf {S} :\delta \mathbf {E} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\delta \mathbf {E} _{\mu }~dV,}$
(9.25)

where ${\textstyle \mathbf {S} }$ and ${\textstyle \mathbf {S} _{\mu }}$ are the macroscopic and microscopic stress tensor, respectively.

Using 9.24, the principle is rewritten as

 ${\displaystyle \mathbf {S} :\delta \mathbf {E} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {S} _{\mu }~dV:\delta \mathbf {E} \left(\mathbf {X} _{o}\right)+{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\delta \mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right)~dV.}$
(9.26)

Taking the macroscopic stress tensor ${\textstyle \mathbf {S} }$ as the volume average of the microstructural stress tensor ${\textstyle \mathbf {S} _{\mu }}$ in the RVE domain, which is similar to the first average relation (see 9.10)

 ${\displaystyle \mathbf {S} \left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\equiv {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {S} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV,}$
(9.27)

equation 9.26 will be satisfied if

 ${\displaystyle \int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\delta \mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right)~dV=\int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\nabla ^{s}\mathbf {\delta w} ~dV=0,}$
(9.28)

Therefore, the RVE's variational equilibrium equation is

 ${\displaystyle \int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\nabla ^{s}\mathbf {\delta w} ~dV=0,}$
(9.29)

which must be satisfied for any kinematically admissible displacement fluctuation field ${\textstyle \mathbf {w} }$ (see Section 9.3.2).

It is possible to observe that because of the symmetry of the stress tensor ${\textstyle \mathbf {S} _{\mu }}$ it can be proved that ${\textstyle \mathbf {S} _{\mu }:\left(\nabla \mathbf {a} \right)=\mathbf {S} _{\mu }:\left(\nabla \mathbf {a} \right)^{T}}$, where ${\textstyle \mathbf {a} }$ is a first order tensor, the 9.28 also can be rewritten as

 ${\displaystyle \int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\delta \mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right)~dV=\int _{\Omega _{\mu }}\mathbf {S} _{\mu }:\nabla \delta \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV=0.}$
(9.30)

### 9.3.5 Microscopic and macroscopic stress tensor

The homogenized stress tensor in the macroscopic level is given by 9.27. The microscopic stress tensor can be obtained as

 ${\displaystyle {\begin{array}{rcl}\displaystyle \mathbf {S} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)&=&\displaystyle \mathbf {C} _{\mu }\left(\mathbf {X} _{\mu }\right):\mathbf {E} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right),\\[3mm]&=&\displaystyle \mathbf {C} _{\mu }\left(\mathbf {X} _{\mu }\right):\mathbf {E} \left(\mathbf {X} _{o}\right)+\mathbf {C} _{\mu }\left(\mathbf {X} _{\mu }\right):\mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right),\end{array}}}$
(9.31)

where ${\textstyle \mathbf {C} _{\mu }}$ is the material constitutive tensor in the RVE. Then, the macro stress tensor is

 ${\displaystyle \mathbf {S} \left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)=\mathbf {\bar {C}} :\mathbf {E} \left(\mathbf {X} _{o}\right)+{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {C} _{\mu }:\mathbf {E} _{\mu }^{w}\left(\mathbf {X} _{\mu }\right)~dV,}$
(9.32)

where,

 ${\displaystyle \mathbf {\bar {C}} \equiv {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {C} _{\mu }~dV}$
(9.33)

is a constitutive tensor which can be considered a microstructural material property.

Equation 9.32 shows that the stress tensor ${\textstyle \mathbf {S} }$ depends of the macroscopic strain tensor ${\textstyle \mathbf {E} }$ and also, of the microscopic strain tensor ${\textstyle \mathbf {E} _{\mu }^{w}}$, which is obtained with the displacement fluctuation field of the RVE. Moreover, the microscopic position vector ${\textstyle \mathbf {X} _{\mu }}$ does not appear explicitly in the microstructural strain tensor expression (see 9.24). Consequently, this variable does not appear in the microstructural stress tensor either. Therefore, the periodic microstructure around the macro point ${\textstyle \mathbf {X} _{o}}$ does not have to be modeled with its exact dimensions. A non dimensional RVE with the internal distribution and volume fractions of the simple materials is enough to obtain the microscopic strain and stress fields. This is one of the principal advantages of this first-order homogenization approach compared to other multiscale high-order approaches.

On the other hand, it can be observed that the kinematically admissible displacement fluctuation option used to satisfy the boundary condition in the RVE problem affects the final macroscopic stress tensor obtained, as occurs in the Taylor model case. This means that if there is a null displacement fluctuation field in the total RVE domain, the stress tensor ${\textstyle \mathbf {S} }$ obtained only depend of the strain tensor ${\textstyle \mathbf {E} }$ and the constitutive tensor ${\textstyle \mathbf {\bar {C}} }$. In other words, the Taylor model condition returns the classical mixing theory results.

## 9.4 Enhanced-first-order homogenization approach

A new enhanced first order homogenization is proposed in the following, in order to enrich the displacement field of the micro model with second order information available in the macro model. The deformed position of a material point in the RVE given by 9.5, used in the first-order approach, can be improved if the second-order term of 9.3 is included. Then, it is possible to propose a new approximation of the current configuration of the RVE as

 ${\displaystyle \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {x} _{\mu }^{c}+\mathbf {F} \left(\mathbf {X} _{o}\right)\cdot \triangle \mathbf {X} _{\mu }+{\frac {1}{2}}\mathbf {G} \left(\mathbf {X} _{o}\right):\triangle \mathbf {X} _{\mu }\otimes \triangle \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right),}$
(9.34)

and setting the coordinate system's origin as defined in 9.6, the proposed deformed position of the RVE is

 ${\displaystyle \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {F} \left(\mathbf {X} _{o}\right)\cdot \mathbf {X} _{\mu }+{\frac {1}{2}}\mathbf {G} \left(\mathbf {X} _{o}\right):\mathbf {X} _{\mu }\otimes \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right).}$
(9.35)

Therefore, the proposed displacement field ${\textstyle \mathbf {u} _{\mu }}$ on the RVE (see 9.8) can be obtained now as

 ${\displaystyle \mathbf {u} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \left[\mathbf {F} \left(\mathbf {X} _{o}\right)-\mathbf {I} \right]\cdot \mathbf {X} _{\mu }+{\frac {1}{2}}\mathbf {G} \left(\mathbf {X} _{o}\right):\mathbf {X} _{\mu }\otimes \mathbf {X} _{\mu }+\mathbf {w} \left(\mathbf {X} _{\mu }\right).}$
(9.36)

Noting that an extra term appears in the proposed microscopic displacement field by including the gradient of the deformation gradient tensor ${\textstyle \mathbf {G} }$ in 9.34. This extra second-order term is a new linking term between the macroscopic and microscopic scales. The proposed displacement field in the RVE is enhanced because it reaches more information from the macro scale.

### 9.4.1 Kinematically admissible displacement fields and boundary conditions in the RVE

The first of the average postulates (see 9.10) is used again to obtain the admissible displacement fields. The microscopic deformation gradient considering the expression given by 9.35 is

 ${\displaystyle \mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)=\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\cong \mathbf {F} \left(\mathbf {X} _{o}\right)+\mathbf {G} \left(\mathbf {X} _{o}\right)\cdot \mathbf {X} _{\mu }+\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right).}$
(9.37)

And, the volume average of this deformation gradient ${\textstyle \mathbf {F} _{\mu }}$ over the RVE is

 ${\displaystyle {\begin{array}{rcl}\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV&=&\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV,\\[5mm]&=&\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)+\mathbf {G} \left(\mathbf {X} _{o}\right)\cdot {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {X} _{\mu }~dV\\[5mm]&+&\displaystyle {\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV.\end{array}}}$
(9.38)

It can be proved that if the RVE geometry in the reference configuration is originally a cube, as shown in Figure 25, and the position of the origin of the coordinate system is defined at the center of the RVE, then the first moment of volume of the RVE is

 ${\displaystyle \int _{\Omega _{\mu }}\mathbf {X} _{\mu }~dV=\mathbf {0} .}$
(9.39)

Therefore, 9.38 can be rewritten as

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {F} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV,}$
(9.40)

or

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)~dV.}$
(9.41)

Equation (9.41) can be rewritten in terms of surface integrals applying the divergence theorem as

 ${\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA-{\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA.}$
(9.42)

Equation 9.42 is exactly the same than the one obtained for the first-order approach (see 9.15) then, the restrictions on the displacement fluctuation field are the same as those shown in 9.16-9.18, and the different displacement fluctuation fields kinematically admissible presented in Section 9.3.2 are still valid for this enhanced-first-order approach.

#### 9.4.1.1 Extra kinematic restrictions and boundary conditions because of G

The next step consists in obtaining the kinematic restrictions result of including the new term ${\textstyle \mathbf {G} }$ in the microscopic displacement field ${\textstyle \mathbf {u} _{\mu }}$. In other words, some extension of the first average theorem needs to be proposed in term of the gradient of the deformation gradient. In the following, a natural extension for the first average theorem is presented. The main drawback of this proposal is that arrives to restrictions on the derivative displacement fluctuation field and therefore, a high-order problem on the RVE must be considered. To avoid this situation, and continue using the classical first-order boundary value problem on the RVE, the alternative extension proposed by Kouznetsova [115] is also presented.

## Natural extension of the first average theorem

The first natural possibility for this extension could be

 ${\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {G} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV.}$
(9.43)

Note that 9.43 is similar to 9.10 but in this case, the volume average of the microstructural gradient of the deformation gradient tensor ${\textstyle \mathbf {G} _{\mu }}$ over the RVE must be equal to the macroscopic gradient of the deformation gradient ${\textstyle \mathbf {G} }$ in the considered point ${\textstyle \mathbf {X} _{o}}$.

Considering 9.37 and 9.4 the gradient of the deformation gradient in the microstructural scale is

 ${\displaystyle \mathbf {G} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)=\nabla \left(\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\right)\cong \mathbf {G} \left(\mathbf {X} _{o}\right)+\nabla \left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\right).}$
(9.44)

Using 9.44 and taking the volume average over the RVE it is possible to obtain

 ${\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\mathbf {G} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\right)~dV,}$
(9.45)

or

 ${\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \left(\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\right)~dV-{\frac {1}{V_{\mu }}}\int _{\Omega _{\mu }}\nabla \left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\right)~dV.}$
(9.46)

And, applying the divergence theorem in the last expression

 ${\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)={\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA-{\frac {1}{V_{\mu }}}\int _{\partial \Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA.}$
(9.47)

Similarly as in the first-order approach, to satisfy the proposed extension of the first average theorem, the integrals that depend of the displacement fluctuation in 9.45 and 9.47 must vanish, then

 ${\displaystyle \int _{\Omega _{\mu }}\nabla \left(\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\right)~dV=\mathbf {0} ,}$
(9.48)

and

 ${\displaystyle \int _{\partial \Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {N} ~dA=\mathbf {0} .}$
(9.49)

The last expression represents a new extra integral restriction on the derivative displacement fluctuation field. Taking the same consideration than before regarding the geometry of the RVE (see Figure 26), the boundary integration in 9.49 can be splitted in

 ${\displaystyle {\begin{array}{l}(\int _{\mathbf {N} _{X}^{+}}\nabla \mathbf {w} ~dA_{yz}-\int _{\mathbf {N} _{X}^{-}}\nabla \mathbf {w} ~dA_{yz})\otimes \mathbf {N} _{X}^{+}\\[5mm]+(\int _{\mathbf {N} _{Y}^{+}}\nabla \mathbf {w} ~dA_{xz}-\int _{\mathbf {N} _{Y}^{-}}\nabla \mathbf {w} ~dA_{xz})\otimes \mathbf {N} _{Y}^{+}\\[5mm]+(\int _{\mathbf {N} _{Z}^{+}}\nabla \mathbf {w} ~dA_{xy}-\int _{\mathbf {N} _{Z}^{-}}\nabla \mathbf {w} ~dA_{xy})\otimes \mathbf {N} _{Z}^{+}=\mathbf {0} .\end{array}}}$
(9.50)
 Figure 27: Normal vectors to the lines in the YZ surface of the Cubic RVE.

Some components of the integrals can also be rewritten in terms of line-boundary integrals applying the divergence theorem. For example, if the first left integral in the first term in 9.50 is taken, the line-boundary of this surface integral can be separated in four different lines, two perpendiculars to Y axis, and the other two perpendiculars to Z axis, as it is shown in Figure 27. Because of the RVE geometry considered, these lines boundary have the property of ${\textstyle \mathbf {N} _{X\mid Y}^{-}=-\mathbf {N} _{X\mid Y}^{+}}$ and ${\textstyle \mathbf {N} _{X\mid Z}^{-}=-\mathbf {N} _{X\mid Z}^{+}}$. Then, with this information the considered integral can be rewritten as

 ${\displaystyle {\begin{array}{l}\int _{\mathbf {N} _{X}^{+}}\nabla \mathbf {w} ~dA_{yz}{\hbox{=}}\int _{\mathbf {N} _{X}^{+}}\nabla _{X}\mathbf {w} ~dA_{yz}\\[5mm]+(\int _{\mathbf {N} _{X\mid Y}^{+}}\mathbf {w} ~dL_{z}-\int _{\mathbf {N} _{X\mid Y}^{-}}\mathbf {w} ~dL_{z})\otimes \mathbf {N} _{X\mid Y}^{+}\\[5mm]+(\int _{\mathbf {N} _{X\mid Z}^{+}}\mathbf {w} ~dL_{y}-\int _{\mathbf {N} _{X\mid Z}^{-}}\mathbf {w} ~dL_{y})\otimes \mathbf {N} _{X\mid Z}^{+},\end{array}}}$
(9.51)

where ${\textstyle \nabla _{X}}$ represents the derivative with respect to the X axis, this term cannot be reduced to a line-integral using the divergence theorem. It can be seen that when the Periodic boundary fluctuations condition is the kinematically admissible option used for the displacement fluctuation field on the RVE, the two right terms on 9.51 are satisfied directly. The points on the opposing lines are a pair boundary points that have the same displacement fluctuation because of the kinematic condition imposed. Applying this same procedure to the rest of the terms of expression 9.50, this equation can be rewritten as

 ${\displaystyle {\begin{array}{l}(\int _{\mathbf {N} _{X}^{+}}\nabla _{X}\mathbf {w} ~dA_{yz}-\int _{\mathbf {N} _{X}^{-}}\nabla _{X}\mathbf {w} ~dA_{yz})\otimes \mathbf {N} _{X}^{+}\\[5mm]+(\int _{\mathbf {N} _{Y}^{+}}\nabla _{Y}\mathbf {w} ~dA_{xz}-\int _{\mathbf {N} _{Y}^{-}}\nabla _{Y}\mathbf {w} ~dA_{xz})\otimes \mathbf {N} _{Y}^{+}\\[5mm]+(\int _{\mathbf {N} _{Z}^{+}}\nabla _{Z}\mathbf {w} ~dA_{xy}-\int _{\mathbf {N} _{Z}^{-}}\nabla _{Z}\mathbf {w} ~dA_{xy})\otimes \mathbf {N} _{Z}^{+}=\mathbf {0} .\end{array}}}$
(9.52)

The previous expression represents an extra restriction on the displacement fluctuation field that makes it kinematically admissible in the RVE. A possible set of boundary conditions that satisfies this restriction is

 ${\displaystyle {\begin{array}{l}{\displaystyle \int _{\mathbf {N} _{X}^{+}}\nabla _{X}\mathbf {w} ~dA_{yz}=\int _{\mathbf {N} _{X}^{-}}\nabla _{X}\mathbf {w} ~dA_{yz}},\\[5mm]{\displaystyle \int _{\mathbf {N} _{Y}^{+}}\nabla _{Y}\mathbf {w} ~dA_{xz}=\int _{\mathbf {N} _{Y}^{-}}\nabla _{Y}\mathbf {w} ~dA_{xz}},\\[5mm]{\displaystyle \int _{\mathbf {N} _{Z}^{+}}\nabla _{Z}\mathbf {w} ~dA_{xy}=\int _{\mathbf {N} _{Z}^{-}}\nabla _{Z}\mathbf {w} ~dA_{xy}}.\end{array}}}$
(9.53)

Equation 9.52 is analogous to 9.18 but it is written in terms of derived displacement fluctuation field, in this case, on the normal direction of the pair surfaces (see Figure 26). Therefore, to satisfy any kinematic restriction, as for example 9.53, obtained from 9.52 a high-order problem on the microscopic scale must be considered because the restriction of the displacement fluctuation field is written on its derivate.

## Alternative extension of the first average theorem

An alternative to the proposed extension of the averaging theorem given by 9.43 should be found to keep a classical boundary value problem on the microstructural RVE problem. With this aim Kouznetsova [115] proposed another alternative extension of the first average theorem. The proposed condition imposes that the second moment of area of the deformed RVE, given in terms of the microscopic displacements, must be equal to the second moment of area of the RVE expressed in terms of macroscopic deformation variables [121]. Considering the above, the expression given by 9.37 is multiplied by ${\textstyle \mathbf {X} _{\mu }}$ and integrated over the RVE volume to obtain

 ${\displaystyle {\begin{array}{rcl}{\displaystyle \int _{\Omega _{\mu }}\nabla _{0}\mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\otimes \mathbf {X} _{\mu }~dV}&=&{\displaystyle \mathbf {F} \left(\mathbf {X} _{o}\right)\otimes \int _{\Omega _{\mu }}\mathbf {X} _{\mu }~dV}\\[4mm]&+&{\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)\cdot \int _{\Omega _{\mu }}\mathbf {X} _{\mu }\otimes \mathbf {X} _{\mu }~dV}\\[4mm]&+&{\displaystyle \int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {X} _{\mu }~dV}.\end{array}}}$
(9.54)

Knowing that the first moment of volume of the undeformed RVE is zero (see 9.39), and defining the second moment of volume of the undeformed RVE as ${\textstyle \mathbf {J} =\int _{\Omega _{\mu }}\mathbf {X} _{\mu }\otimes \mathbf {X} _{\mu }~dV}$. Equation 9.54 can be rewritten as

 ${\displaystyle \mathbf {G} \left(\mathbf {X} _{o}\right)\cdot \mathbf {J} =\int _{\Omega _{\mu }}\nabla \mathbf {x} _{\mu }\left(\mathbf {X} _{o},\mathbf {X} _{\mu }\right)\otimes \mathbf {X} _{\mu }~dV-\int _{\Omega _{\mu }}\nabla \mathbf {w} \left(\mathbf {X} _{\mu }\right)\otimes \mathbf {X} _{\mu }~dV,}$
(9.55)

replacing the following relationships