Numerical homogenization of piezoelectric textiles with electrospun fibers for energy harvesting

Piezoelectric effects are exploited in an increasing number of microand nano-electro-mechanical systems. In particular, energy harvesting devices convert ambient energy (i.e. mechanical pressure) into electrical energy and their study is nowadays a very important and challenging field of research. In this paper, the attention is focused on piezoelectric textiles. Due to the importance of computational modeling to understand the influence that micro-scale geometry and constitutive variables have on the macroscopic behavior, a homogenization strategy is developed. The macroscopic structure behaviour is obtained defining a reference volume element (RVE) at the micro-scale. The geometry of the RVE is based on the microstructural properties of the material under consideration and consists in piezoelectric polymeric nano-fibers subjected to electromechanical contact constraints. This paper outlines theory and numerical implementation issues for the homogenization procedure. Moreover, within this approach the average response resulting from the analysis of different fiber configurations at the microscale is determined providing a multiphysics constitutive model for the macro-scale.


INTRODUCTION
lectrospinning is a simple and versatile method for generating ultrathin fibers from a rich variety of materials that include polymers, composites and ceramics [1].This non-mechanical, electrostatic technique involves the use of a high voltage electrostatic field to charge the surface of a polymer solution droplet and thus to induce the ejection of a liquid jet through a spinneret.Nowadays, nanofiber technology is opening up new scenarios in several industrial fields.Applications range from energy harvesting technologies [2,3] to tissue and biomedical engineering, aerospace materials, and device integration with architectural or design components.Nanofibers can increase the performance of traditional materials and allow engineering optimization of existing textiles and fabrics and even development of new E materials.From a commercial perspective the term "nano" describes a diameter of the fiber below one micron.However, the most important properties of these materials tend to manifest at a scale below 500 nanometers.Materials characterized by nanofibers at the microlevel show high surface area and superior mechanical, electric, magnetic properties [4].In this framework, polyvinylidene fluoride (PVDF) is a fluoropolymer known for its strength, chemical resistance, thermal resistance and piezoelectric properties.In particular PVDF nanofibers are suitable for numerous applications such as filtration, coating, sensors, and energy generators.From the chemical point of view, the PVDF material is built joining chains of CH2CF2, where C indicates the carbon, H the hydrogen and F the fluorine atoms.It is produced in large thin clear sheets and through a stretching and poling process it is possible to give piezoelectric properties to the resulting thin layer.The stretch direction is the direction along the sheet in which most of the carbon chains run.The hydrogen atoms, which have a net positive charge, and the fluorine atoms, which have a net negative charge, end up on opposite sides of the sheet.This creates a pole direction that is either oriented to the top or bottom of the sheet.When an electric field E is applied across the sheets, they either contract in thickness and expand along the stretch direction or expand in thickness and contract along the stretch direction depending on which way the field is applied.This is due to the physical nature of the positive hydrogen atoms attracted by the negative side of the electric field and repelled by the positive side of the electric field.Scanning electron micrographs of PVDF layer architectures at the microscale show that, depending on the process parameters, the resulting microstructure can range from a fully random geometry to excellent mutual alignment of fibers [5].Several analytical and computational multiscale approaches have been developed in order to predict the macroscale properties of heterogeneous materials at the lower scale(s), [6][7][8][9][10][11].Although most of these efforts have been devoted to continuum mechanics [12,13], some applications to multiphysics problems are also available [14,15].A few of these concern electromechanically coupled problems such as in the case of piezoelectricity [16].Moreover, although several macroscale formulations were developed for piezoelectric shell elements [17][18][19], computational homogenization of shells is only recently receiving major attention [20,21].To the best of our knowledge, such approaches have not yet been proposed for piezoelectric shells.In this framework, the objective of this paper is to model the behaviour of the aforementioned PVDF sheet by a multiscale and multiphysics approach.This requires the definition of a representative volume element (RVE) at the microscale, the formulation and solution of a microscale boundary value problem (BVP), and the development of a suitable micro-macro scale transition.In the first part of the paper a microscale RVE element is defined and some details about the formulation of suitable electromechanical contact laws to describe the interaction among the fibers are provided.The second part describes the kinematic behaviour of a piezoelectric shell and the multiscale approach.Finally, based on the presented framework, several RVE geometries are analysed and the homogenized coefficients determined.

MICROSCALE FORMULATION AND HOMOGENIZATION PROCEDURE
he microscopic length scale of the RVE (micrometers) is several orders of magnitude smaller than the macroscopic dimensions (centimeters).Hence, two models are developed: one at the micro-level and one at the macro-level.The problem is how to couple these models and which boundary conditions to apply to the micromodel.At the microscale, the RVE consists in piezoelectric polymer fibers that feature a linear piezoelastic constitutive behavior and are subjected to electromechanical contact constraints.The governing equations are the Navier equations and the strain-displacement relations for the mechanical field and Gauss and Faraday laws for the electrostatic field.Moreover, the constitutive equations read: where ijkl C , ikl e , and  ik  are respectively the elastic, piezoelectric, and permittivity constants, whereas , ij ij S T are the strain and stress components and , i i D E are the electric displacement and the electric field components, respectively.The interaction among the fibers is described defining a 3D electromechanical frictionless contact law and implementing an element with the following main characteristics:  the contact formulation is based on the master-slave concept;  Bézier patches are used for smoothing of the master surface;  the impenetrability condition is extended to the electromechanical setting by imposing equality of the electric T potential in case of closed contact.The electromechanical constraints are regularized with the penalty method.For each slave node, the normal gap is computed as: where s x is the position vector of the slave node, m x is the position vector of its normal (i.e.minimum distance) projection point onto the master surface, and n is the outer normal to the master surface at the projection point.The sign of the measured gap is used to discriminate between active and inactive contact conditions, a negative value of the gap leading to active contact.The electric field requires the definition of the contact electric potential jump: where  s and  m are the electric potential values in the slave node and in its projection point on the master surface.A tensor product representation of one-dimensional Bézier polynomials is used to interpolate the master surfaces in the contact interface for a three-dimensional problem according to the relation: where  kl is the potential evaluated at 16 control points as a function of the potential at the auxiliary points kj and at the master nodes  m ij , see Fig. 1.According to standard finite element techniques, the global energy of the system is obtained by adding to the variation of the energy potential representing the continuum behaviour the virtual work associated to the electromechanical contact contribution provided by the active contact elements.Performing the first and second variation of the global energy the global set of equations is obtained.
When an external load is applied to the RVE, the stress, strain and electric fields in the microstructure will show large gradients due to the microstructural heterogeneity.However, due to the differences in scale, the microstructural electric/deformation field around a macroscopic point will be approximately the same as the electric/deformation field around neighbouring points.The repetitive deformations justify the assumption of local periodicity, meaning that the microstructure can be thought as repeating itself near a macroscopic point.However, the microstructure itself may differ from one macroscopic point to another.The repetitive microstructural generalized deformations suggest that macroscopic stresses and strains around a certain macroscopic point can be found by averaging microstructural stresses and strains in a small representative area of the microstructure attributed to that point.This allows finding a globally homogeneous medium equivalent to the original composite, where the equivalence is intended in an energetic sense as per Hill's balance condition [12,15,16].Formulated for the electromechanical problem, Hill's criterion in differential form reads: and requires that the macroscopic volume average of the variation of work performed on the RVE is equal to the local variation of work on the macroscale.In the previous equation: ij T , ij S , i D and i E represent respectively the average values of stress, strain, electric displacement and electric field components and V is the RVE volume.Hill's lemma leads to the following equations: Fig. 2 shows a two-dimensional RVE.Periodic boundary conditions imply that, on two opposite edges, the displacement and the electric potential are equal, the stress vector and the electric displacement vector are opposite.In the implementation, these constraints are enforced by prescribing the primal variables at the corner nodes and using the Lagrange multiplier method.X X ) boundaries this relation is valid in the reference configuration: where p X , p = 1, 2, 4 are the position vectors of the corner nodes 1, 2, and 4 in the undeformed state.Then in the deformed configuration the previous relations lead to where M F is the deformation gradient.Now if the position vectors of the corner nodes in the deformed state are prescribed according to: then the periodic boundary conditions may be rewritten in terms of displacement and electric potential as: where p where the macroscopic (overall, effective) mechanical moduli C , piezoelectric moduli e , and dielectric moduli   , are introduced.The overall material behaviour resulting from the homogenization procedure is nonlinear due to the electromechanical frictionless contact conditions at the microscale.Therefore, the coefficients in Eq. ( 12) are to be considered secant values determined for a given load increment.For the next steps we define the constitutive matrix macro solid

D
and the generalized Green Lagrange strain vector g E of the homogenized solid as: Within this approach only constitutive equations at the RVE scale are required.

MACROSCALE FORMULATION
ased on the analysis of the shell kinematics, see Fig. 3, the Green-Lagrange strains and the electric field components in convective coordinates can be arranged in a generalized strain column vector:

B
In eq. 14 the membrane strain components   and the change of curvature   read: where comma indicates partial derivation, Greek indices take the values 1, 2; 0 , ψ ψ are respectively the current and initial position vectors of the shell middle surface, g is the initial shell director and 0 h is the initial shell thickness.
Moreover we have introduced the quantity: where 1 2 ,   are the natural coordinates of the shell middle surface,   is the thickness stretch, h is the current shell thickness and a is the current normal.Furthermore in eq.14, the shear strain components  take the form: and the electric components  E read: where  is the electric potential.Finally 0 1 33 33 ,   are the constant and linear components of the thickness strain and , E E represent the constant and linear parts of the electric field along the thickness direction.We now introduce the transformation matrix A between the generalized Green Lagrange strain vector of the solid g E and the generalized strain column vector of the shell s E such as  g s E AE with: 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 With some algebra and after integration on the shell thickness the final constitutive equation of the homogenized shell can be recast in the following form:  , , , ,  is the natural coordinates of the shell thickness.The above formulation is valid for the general case of finite strain shell problems.However, as a first step we implemented here only its linearized version, whereas the finite strain implementation is the subject of ongoing research.

RESULTS
s follows, the presented computational procedure is applied to simple test problems.For the general case of a solid RVE with random fibers distribution, see Fig. 4, the final material constitutive equations macro solid

D
will be those of an anisotropic piezoelectric solid, see eq.24. (24) For the evaluation of the effective properties, suitable boundary conditions have to be applied to the unit cell in such a way that, apart from one component of the strain/electric field vector, all other components are equal to zero [14,15].Then each effective coefficient can be easily determined by multiplying the corresponding row of the material matrix by the strain/electric field vector.Once macro solid D is obtained, coefficients of macro shell D can be easily derived using eq.22.In Fig. 4 the mesh used during the analyses are illustrated.Each RVE is assumed to have a dimension A x B with A=B=10  m and height 4  m.In Fig. 10 (a-e) the coefficients as obtained from the homogenization procedure (black dot) are compared with the corresponding coefficients derived considering a bulk RVE with the properties of the fibers (red dot) and the coefficients obtained scaling the bulk values to take into account the volume of voids in the RVE.Both the fibers arrangement and the void fraction lead to an homogenized material behaviour more flexible than the bulk.This means that the final material will present improved capabilities when used to build sensors and nanogenerators since a lower force will cause a higher deformation in the material resulting in an increased output voltage.

CONCLUSIONS
his paper proposes a general formulation for homogenization of the electromechanical behaviour of piezoelectric textile nanogenerators built assembling PVDF nanofibers using an electrospinning process.In particular the effects of microstructure geometry and fiber distribution in the RVE are investigated.Despite the resulting fibrous material has a main polarization along the longitudinal axis of the fibers, the interactions among fibers lead to a complex three-dimensional distribution of stress, strain and electric potential in the RVE.The results shed light on the homogenized response of piezoelectric textiles where the resulting material can be described with an anisotropic piezoelectric constitutive matrix with further coupling coefficients that were equal to zero at the micro level (fiber scale).Based on the proposed approach, it is possible to calculate anisotropic material constants of an equivalent homogenized piezoelectric solid and shell.Some presented results demonstrate the capability of the developed procedure and algorithms.Further studies are needed for a full characterization of the macroscopic material behaviour aiming both at introducing in the numerical model more reliable electromechanical contact laws based on experimental results under development and at implementing a full coupling between the micro and macro scales in the framework of FE 2 methods.
are the Bernstein polynomials and kl d the coordinates of the control points[22].With the same procedure, to interpolate the electric potential on the master surface the following relation is introduced:

Figure 1 :
Figure 1: Electromechanical contact elements: a) General master-slave concept b) Node to surface discretization c) Smoothing with Bezier patches

Figure 2 :
Figure 2: Enforcing periodic boundary conditions During the initial periodicity of the RVE, for every respective pair of nodes on the top-bottom ( 4 3 1 2 ,   X X ) and right-

u
and  p , p = 1, 2, 4 are the position vectors and the electric potential of the corner nodes 1, 2, and 4 in the undeformed state and ( position vectors and electric potential for every respective pair of nodes on the top-bottom and right-left boundaries of the RVE.Once the RVE problem is solved, the macroscopic material properties are determined from the homogenization procedure and the final constitutive equations read:

n
are the membrane forces,  m are the bending moments,  q are the shear forces, i d the constant and linear components of membrane force and electric displacement in the thickness direction.Moreover 3

Figure 4 : 11 
Figure 4: RVE mesh and boundary conditions.Each fiber is considered as a linear piezoelastic solid and is discretized with linear 8-node brick elements.The fibers are polarized in the direction parallel to the fiber longitudinal axis, i.e. axis 3 in the next figures.Frictionless electromechanical

Figure 5 :
Figure 5: Contours of displacement, stress, electric field and potential: traction in direction 1.

Figure 6 :
Figure 6: Contours of displacement, stress, electric field and potential: compression in direction 1

Figure 7 :
Figure 7: Contours of displacement, stress, electric field and potential: compression in direction 1

Figure 8 :
Figure 8: Contours of displacement, stress, electric field and potential: compression in direction 2

Figure 9 :Figure 10 :
Figure 9: Contours of displacement, stress, electric field and potential: compression in direction 2