Microsoft Word - numero_29_art_6 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 49 Focussed on: Computational Mechanics and Mechanics of Materials in Italy Numerical homogenization of piezoelectric textiles with electrospun fibers for energy harvesting C. Maruccio Università del Salento claudio.maruccio@unisalento.it L. De Lorenzis Technische Universität Braunschweig l.delorenzis@tu-braunschweig.de ABSTRACT. Piezoelectric effects are exploited in an increasing number of micro- and 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. KEYWORDS. Electromechanical Coupling; Multiphysics Modeling; Multiscale Modeling; Shell elements; Energy harvesting. 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 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 50 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-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-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 micro- model. 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: ) )      ij ijkl kl kij k i ikl kl ik k a T C S e E b D e S E (1) where ijklC , ikle , and ik are respectively the elastic, piezoelectric, and permittivity constants, whereas , ij ijS T are the strain and stress components and , i iD 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 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 51 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: ( )  s m Ng x x n (2) where sx is the position vector of the slave node, mx 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: ( )   s mg (3) 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:  1 2 1 2 0 0 , ( ) ( )        m m m m kl k l k l x B Bd (4) where 1 ( )m kB and  2 m lB are the Bernstein polynomials and kld 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:  1 2 1 2 0 0 , ( ) ( )          m m m m kl k l k l B B (5) 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 mij , see Fig.1. a) b) c) Figure 1: Electromechanical contact elements: a) General master-slave concept b) Node to surface discretization c) Smoothing with Bezier patches 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 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 52 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: 1 1       ij ij i i ij ij i iT S D E T S dV D E dV V V (6) 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: ijT , ijS , iD and iE 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: 1 1 1 1         ij ij ij ij i i i i T T dV S S dV V V D DdV E E dV V V (7) 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. 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– left ( 1 4 2 3, X X ) boundaries this relation is valid in the reference configuration: 4 3 1 2 4 1   X X X X 1 4 2 3 2 1   X X X X (8) where pX , 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  4 3 1 2 4 1   Mx x F X X  1 4 2 3 2 1   Mx x F X X (9) where MF is the deformation gradient. Now if the position vectors of the corner nodes in the deformed state are prescribed according to: p M px F X with p = 1, 2, 4 (10) then the periodic boundary conditions may be rewritten in terms of displacement and electric potential as: C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 53         4 3 1 2 4 1 4 3 1 2 4 1 1 4 2 3 2 1 1 4 2 3 2 1                             u u u u u u u u (11) where pu 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 ( 4 3 1 2 1 4 2 3, , ,   u u u u ), ( 4 3 1 2 1 4 2 3, , ,       ) are 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:                 TC eT S D Ee  (12) 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 solidD and the generalized Green Lagrange strain vector gE of the homogenized solid as:         macro solid TC e D e         g S E E (13) 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: 0 1 0 1 11 22 12 11 22 12 1 2 1 2 33 33 3 3, , 2 , , , 2 , , , , , , , ,             T s E E E EE (14) Figure 3: Shell kinematics. B C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 54 In eq. 14 the membrane strain components  and the change of curvature  read: , , 0, 0, 1 ( ) 2        ψ ψ ψ ψ (15) 0 0 , , , , 0, , 0, , 1 ( ) 2 2 2                 h h ψ d ψ d ψ g ψ g (16) 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 0h is the initial shell thickness. Moreover we have introduced the quantity: 0 1 2 1 2( , ) ( , ) 2      h d a (17) where 1 2,  are the natural coordinates of the shell middle surface,  1 2 0 ,    h h 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: 0 , 0,2           h ψ d ψ g (18) and the electric components E read:        E (19) where  is the electric potential. Finally 0 1 33 33,   are the constant and linear components of the thickness strain and 0 1 3 3, 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 gE and the generalized strain column vector of the shell sE such as g s E AE with: 3 3 3 3 3 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                                  A (20) 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:  macro shell sL D E (21) with C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 55 1 2 3 1 2    macro T macro shell solidh dD A D A (22) and 0 1 0 1 11 22 12 11 22 12 1 2 1 2 33 33 3 3, , , , , , , , , , , , ,       T n n n m m m q q d d n n d dL (23) where n are the membrane forces, m are the bending moments, q are the shear forces, id are the electric displacement and 0 1 33 33,n n , 0 1 3 3, d d are the constant and linear components of membrane force and electric displacement in the thickness direction. Moreover 3 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 solidD 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 solidD is obtained, coefficients of macro shellD can be easily derived using eq. 22. In Fig.4 the mesh used during the analyses are illustrated. (a) (b) (c) (d) 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 (d) (a) A B (c) (b) A C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 56 contact constraints are enforced at the interface between the fibers using the electromechanical contact formulation discussed in section 2. Figures from 5 to 9 provide the contours of displacements, stresses electric field and potential for special cases of boundary conditions applied to the RVE as a function of the fiber arrangement in the material at the microscale. Moreover Fig. 10 provides as an example a comparison between some of the resulting homogenized coefficients as a function of the RVE geometry, namely 11C , 12C , 13C , 31e . In particular for each fiber we assumed for the elastic modulus E=2.5 GPa and the Poisson ratio 0.3  , while the piezoelectric strain coefficients d31, d32, d33 are 20x10-12, 3x10-12, 33 x10-12 m/V (leading to piezoelectric stress coefficients e31, e32, e33 equal to 0.03077, -0.0057, -0.075 C/m2) and the permittivity coefficients 11 , 22 , 33 are equal to 12 0 , 0 being the vacuum permittivity. 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. Figure 5: Contours of displacement, stress, electric field and potential: traction in direction 1. Figure 6: Contours of displacement, stress, electric field and potential: compression in direction 1 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 57 Figure 7: Contours of displacement, stress, electric field and potential: compression in direction 1 Figure 8: Contours of displacement, stress, electric field and potential: compression in direction 2 Figure 9: Contours of displacement, stress, electric field and potential: compression in direction 2 C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 58 a) b) c) d) e) Figure 10: Homogenized RVE coefficients. C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 59 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 FE2 methods. ACKNOWLEDGEMENTS he authors acknowledge the support from the Italian MIUR through the project FIRB Futuro in Ricerca 2010 "Structural mechanics models for renewable energy applications" (RBFR107AKG) and from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013), ERC Starting Grant INTERFACES (grant agreement n. 279439). REFERENCES [1] Dzenis, Y., Spinning continuous fibers for nanotechnology, Science 304 (5679) (2004) 1917–1919. [2] Gao, Y., Wang, Z.L., Electrostatic potential in a bent piezoelectric nanowire, The fundamental theory of nanogenerator and nanopiezotronics, Nano Lett, 7 (2007) 2499-2505. [3] Persano, L., Dagdeviren, C., Su, Y., Zhang, Y., Girardo, S., Pisignano, D., Huang, Y., Rogers, J. A., High performance piezoelectric devices based on aligned arrays of nanofibers of PVDF, Nature Communications, 1633 (2013) 4. [4] Pisignano, D., Polymer nanofibers. Cambridge: Royal Society of Chemistry, (2013). [5] Maruccio, C., De Lorenzis, L., Persano, L., Pisignano, D., Computational homogenization of fibrous piezoelectric materials. Submitted. [6] Huet, C., Application of variational concepts to size effects in elastic heterogeneous bodies, J. Mech. Phys. Solids, 38(6) (1990) 813–841. [7] Terada, K., Hori, M., Kyoya T., Kikuchi N., Simulation of the multi-scale convergence in computational homogenization approach, Int. J. Solids Struct., 37 (2000) 2285–2311. [8] Kanit, T., Forest, S., Galliet, I., Mounoury V., Jeulin, D., Determination of the size of the representative volume element for random composites: Statistical and numerical approach, Int. J. Solids Struct., 40 (2003) 3647–3679. [9] Bisegna, P., Luciano, R., Bounds on the overall properties of composites with debonded frictionless interfaces. Mechanics of materials, 28 (1998) 23-32. [10] Kouznetsova, V., Brekelmans, W. A. M., Baaijens F. P. T., An approach to micro-macro modeling of heterogeneous materials. Comput. Mech., 27 (2001) 37-48. [11] Geers, M. G. D., Kouznetsova, V. G., Brekelmans, W. A. M., Multi-scale first-order and second-order computational homogenisation of microstructures towards continua, Int. J. Multiscale Comput. Eng., 1 (2003) 371-386. [12] Kouznetsova, V., Geers, M. G. D., Brekelmans, W. A. M., Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenisation scheme, Int. J. Numer. Meth. Eng., 54 (2002) 1235-1260. [13] Miehe, C., Computational micro-to-macro transitions for discretized microstructures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy, Comput. Methods Appl. Mech. Eng., 192 (2003) 559–91. T T C. Maruccio et alii, Frattura ed Integrità Strutturale, 29 (2014) 49-60; DOI: 10.3221/IGF-ESIS.29.06 60 [14] Berger, H., Gabbert, U., Koeppe, H., Rodriguez-Ramos, R., Bravo-Castillero, J., Diaz, G. R., Otero, J. A., Maugin, G. A., Finite element and asymptotic homogenization methods applied to smart composite materials, Comput. Mech., 33 (2003) 61-7. [15] Berger, H., Kari, S., Gabbert, U., Rodriguez-Ramos, R., Guinovart-Diaz, R., Otero, J. A., Bravo-Castillero J., An analytical and numerical approach for calculating effective material coefficients of piezoelectric fiber composites, Int. J. Solids Struct., 42 (2004) 5692-714. [16] Schroeder J., Keip, M., Two-scale homogenization of electromechanically coupled boundary value problems - Consistent linearization and applications, Computational Mechanics, 50(2) (2012) 229-244. [17] Schulz, K., Klinkel, S., Wagner, W., A finite element formulation for piezoelectric shell structures considering geometrical and material non-linearities, Int. J. Numerical Methods in Engineering, 87 (2011) 491–520. [18] Klinkel, S., Wagner, W., A piezoelectric solid shell element based on a mixed variational formulation for geometrically linear and nonlinear applications. Computers and Structures, 86 (2008) 38–46. [19] Klinkel, S., Gruttmann, F., Wagner,W., A mixed shell formulation accounting for thickness strains and finite strain 3d material models, Int. J. Numerical Methods in Engineering, 75 (2008) 945–970. [20] Fillep, S., Mergheim, J., Steinmann, P., Computational modelling and homogenization of technical textiles. Eng. Structures, 50 (2013) 68–73. [21] Coenen, E. W., Kouznetsova, V. G., Geers, M. G. D., Computational homogenization for heterogeneous thin sheets. International Journal for Numerical Methods in Engineering, 83 (2010) 1180-1205. [22] Lengiewicz, J., Korelc, J., Stupkiewicz, S., Automation of finite element formulations for large deformation contact problems. International Journal for Numerical Methods in Engineering, 85 (2011) 1252-1279.