Computational Homogenization of Fibrous Piezoelectric Materials

Claudio Maruccio, Laura De Lorenzis, Luana Persano, Dario Pisignano

Dipartimento di Ingegneria dell’Innovazione, Universit del Salento, Via Monteroni, Lecce, Italy. E-mail:

Institut f r Angewandte Mechanik, Technische Universit t Braunschweig, Bienroder Weg, Braunschweig, Germany. E-mail:

National Nanotechnology Laboratory, CNR-Istituto Nanoscienze, via Arnesano, Lecce, Italy. E-mail:

Dipartimento di Matematica e Fisica “E. De Giorgi”, Universit del Salento, via Monteroni, Lecce, Italy. E-mail:


Computational homogenization, electromechanical contact, multiphysics modeling, multiscale modeling, polymer nanofibers, nonlinear piezoelectricity.


Flexible piezoelectric devices made of polymeric materials are widely used for micro- and nano-electro-mechanical systems. In particular, numerous recent applications concern energy harvesting. Due to the importance of computational modeling to understand the influence that microscale geometry and constitutive variables exert on the macroscopic behavior, a numerical approach is developed here for multiscale and multiphysics modeling of thin piezoelectric sheets made of aligned arrays of polymeric nanofibers, manufactured by electrospinning. At the microscale, the representative volume element consists in piezoelectric polymeric nanofibers, assumed to feature a piezoelastic behavior and subjected to electromechanical contact constraints. The latter are incorporated into the virtual work equations by formulating suitable electric, mechanical and coupling potentials and the constraints are enforced by using the penalty method. From the solution of the micro-scale boundary value problem, a suitable scale transition procedure leads to identifying the performance of a macroscopic thin piezoelectric shell element.

1 Introduction

The discovery of appreciable piezoelectricity on ceramic and polymeric materials such as ZnO and Polyvinylidene fluoride (PVDF) has led to the development of a series of nanowire-based piezoelectric nanogenerators [10, 4, 5]. These piezoelectric devices are attractive for several technological applications, most notably for mechanical energy harvesting and pressure/force sensors. The impact of these technologies e.g. for powering small electronic devices and, in perspective, coupling with sensing and photonic platforms, is significant and foreseeably going to grow. In particular, PVDF is a polymeric piezoelectric material with good piezoelectric and mechanical properties. For comparison, the most important piezoelectric coefficients are for ZnO: =12 pC/N, =-4.7 pC/N and =-12 pC/N, and for bulk PVDF =-30 pC/N, =23 pC/N, =0 pC/N. Recent experimental results showed that the piezoelectric coefficient of PVDF in the form of nanofibers is around -40 pC/N, i.e. is higher in absolute value than that of ZnO and bulk PVDF. In particular, considering PVDF in the form of nanofiber arrays an increase up to 20-30 is found as a function of the alignment of the fibers at the microscale [19, 20, 35]. Moreover, the high flexibility typical of polymers allows very high strains to be applied to a PVDF sheet or filament, thus high piezoelectric potentials can be expected. Electrospinning technologies, which exploit the elongation of electrified jets to produce polymer nanofibers, are attracting an increasing attention in this framework, because they are especially effective, cheap and allow good amounts of functional fibers to be realized in continuous runs [24, 18, 29, 23]. In particular, a few electrospinning methods were applied to piezoelectric polymers to realize active, eventually self-poled piezoresponsive fibers. For instance, far-field electrospinning was recently used [22, 35] to produce piezoelectric thin sheets of PVDF nanofibers aligned in uniaxial arrays. This approach allows macroscopic samples (sheets) of aligned nanofibers to be obtained, as well as different microscale fiber architectures to be produced which may vary in a controlled way from a completely randomly oriented to a fully aligned ensemble depending on the process parameters and specific experimental set-up, as shown in Fig. 1. During the fabrication process we observed that welding takes place between adjacent fibers which leads the assembly of nanofibers to behave as a structure. The stretching force and a strong electric field applied between the nozzle and the collecting substrate surface pole the PVDF nanofibers into the phase. The resulting PVDF sheet can then be directly applied between two electrodes on a flexible plastic substrate for piezoelectric output measurement [22, 35]. As shown in Fig. 1e, the material which macroscopically appears in the form of thin sheets is actually made of aligned arrays of polymeric nanofibers. This becomes evident when looking at it one scale below (termed “microscale” in the following).

Due to the importance of understanding the influence that microscale geometry and constitutive variables exert on the macroscopic behavior, the objective of this work is to develop a multiphysics multiscale procedure to obtain the macroscopic piezoelastic behavior of the thin sheets based on the microscale features of the constituent material. In order to predict the macroscale properties of materials and devices featuring heterogeneous properties at the lower scale(s), several analytical and computational multiscale approaches have been developed in the past years [16, 11] to overcome the prohibitive computational expense required for an explicit description of the lower-scale features. Although most of these efforts have been devoted to continuum mechanics [21, 28], some applications to multiphysics problems are also available. A few of these focus specifically on electromechanically coupled problems such as in the case of piezoelectricity [26], and are based either on analytical approaches [3, 25] or on finite element analyses [1, 2]. Although several macroscale formulations were developed for piezoelectric shell elements [27, 12, 13], computational homogenization of shells is only recently receiving major attention [6]. Existing approaches model the representative volume element (RVE) as an ordinary 3D continuum, described by the standard kinematic, equilibrium and constitutive equations. In-plane homogenization is combined with through-thickness integration to obtain the macroscopic generalized stress and moment resultants [9]. To the best of our knowledge, such approaches have not yet been proposed for piezoelectric shells.

In this paper, we present a computational homogenization procedure which derives macroscopic non-linear constitutive laws for piezoelectric shells made of aligned arrays of polymeric nanofibers from the detailed description of the microscale. The motivation is to describe the behavior of the PVDF sheets mentioned earlier and, in a subsequent phase, to optimize their macroscopic piezoelectric response by tailoring their microscale features. The concept is schematically illustrated in Fig. 2.

This paper is organized as follows: Section 2 describes the kinematic behavior of a piezoelectric shell, where displacements and electric potential are the independent fields. In Section 3 a microscale RVE element is defined and the kinematically nonlinear theory of piezoelasticity is briefly reviewed along with its finite element formulation at the microscale. The implementation of suitable frictionless electromechanical contact elements using smoothing techniques is described and details regarding the linearization of the finite element equations at the micro scale are also provided. Section 4 describes the transition between the micro- and macro scales. Finally, in Section 5 the RVE geometry is analyzed to determine the effective material properties of the macroscopic shell. Advanced symbolic computational tools available in the AceGen/AceFEM finite element environment within Mathematica [15, 14] are used throughout this work, with the advantage that the tasks related to the finite element implementation, including linearization of the non-linear governing equations, are largely automated.

2 Macroscale: shell kinematics

In this section we describe the kinematics assumed for the macroscale piezoelectric shell, Fig. 3, following [27], which in turn is largely based on Naghdi’s theory for the mechanical part [33, 36]. For more details, see the original paper [27]. The Green-Lagrange strains and the electric field are derived in convective coordinates. The parameter is defined as the thickness parametric coordinate and with are the in-plane parametric coordinates of the shell middle surface. Thus , with as the domain of the shell middle surface parameterization. We further introduce the convention that indices in Greek letters take the values whereas indices in Latin letters take the values . Moreover, partial derivatives are denoted as follows


2.1 Mechanical field

The position vector in the reference (initial) shell configuration is defined as :


where is the reference position vector of the shell middle surface, is the initial shell thickness, is the initial unit normal, and , with and .

The position vector in the current shell configuration is defined as :


where is the current position vector of the shell middle surface, is the current shell thickness, is the current unit normal, and , with and . In turn, can be expressed as


with as the shell middle surface displacement vector. By introducing the thickness stretch


the current shell director can be rewritten as


Let us now introduce the covariant basis vectors in the reference and current configurations, respectively and , as follows, see Fig. 3




From the covariant bases, the contravariant bases and are induced with




The deformation gradient can be expressed as a function of the basis vectors in the reference and current configurations:


A suitable strain measure is the Green-Lagrange strain tensor, defined as


where is the second-order unit tensor. Substituting eq. (12) into eq. (13) and recalling that the identity tensor is identical to the metric tensor yields


Comparison of eqs. (13) and (14) yields


If the expressions of and in eq. (15) are expanded using eqs. (11), (7) and (8), the in-plane strain components can be expressed as follows:




representing the membrane strain components and


representing the change of curvature. In eq. (16) an additional term quadratic in has been neglected, a usual assumption which delivers accurate results for thin shells such as those considered herein [33, 34]. The shear strain components assume the form




again neglecting an additional term, linear in . Finally, the thickness strain is obtained as


2.2 Electric field

The electric field is given by the gradient of the electric potential , hence in the material configuration its expression is given by


If the piezoelectric material is assumed to be poled in the shell thickness direction, only the difference of electric potential in this direction must be considered, otherwise also the other contributions are required. A common assumption in piezoelectric models is that the electric field is constant through the thickness inside the actuator or sensor. This is in bending dominated situations not correct [27]. Herein, a linear approximation is adopted, which is sufficient to pass the out-of-plane bending patch test [27].

2.3 Generalized strain vector

The Green-Lagrange strain and the electric field components of the piezoelectric solid can be arranged in a generalized strain column vector :


whereas the strain and electric field components of the piezoelectric shell can be arranged in the following vector


where are the constant and linear components of the thickness strain, while represent the constant and linear parts of the electric field in the thickness direction [27]. The relation between the Green-Lagrange strains and the independent shell strains can be stated in compact form by introducing the matrix A such as:




3 The microscale boundary value problem

The material within the RVE consists in piezoelectric polymer fibers with a tight packing arrangement (Fig. 2). The perfectly aligned configuration considered herein is a simplification of the actual geometry, which will be tackled in future extensions. The fibers are partially bonded to each other (as will be described in more detail in Section 4.2) and partially subjected to unilateral electromechanical contact constraints. This section overviews the kinematically nonlinear electroelasticity theory and the contact formulation, whereas the boundary conditions will be discussed in Section 4.2 as they are integral to the scale transition procedure.

3.1 Nonlinear electro-elasticity theory

The starting point for the variational setting of geometrically nonlinear, quasi-static electro-elasticity is the definition of an energy density per unit reference volume [40]:


where is the deformation gradient and is the electric field in the reference configuration given by where is the gradient operator in the reference configuration and is the electric potential. All quantities refer to the microscale as indicated by the superscript m, also used in the following for the same purpose. Based on the above energy density the first Piola Kirchoff stress and the reference dielectric displacement are defined as:


The mechanical equilibrium conditions and the Gauss law of electricity are (in the absence of body forces and charges):


where is the divergence operator. The fourth-order reference elasticity tensor follows as the derivative of the first total Piola Kirchoff stress with respect to the deformation gradient or as part of the hessian of the total energy density :


Likewise the second-order reference dielectricity tensor follows as the derivative of the reference dielectric displacement with respect to the reference electric field or as part of the hessian of the total energy density :


Finally the third-order reference piezoelectricity tensor follows either as the derivative of the reference dielectric displacement with respect to the deformation gradient or alternatively as the derivative of the first Piola Kirchoff stress with respect to the reference electric field:


Objectivity requires that (with a slight notation abuse):


where is the right Cauchy Green tensor, or equivalently:


where is the Green Lagrange strain tensor. Moreover instead of it is often more convenient to use the second Piola Kirchhoff stress . So we can now introduce the fourth-order reference elasticity tensor and the third-order reference piezoelectricity tensor referred to the new variables and as


All constants terms can be grouped in a generalized microscale constitutive matrix :


where and are matrices obtained respectively from the fourth and third order tensors in eqs. (35) and (36) using Voigt notation. Further assuming linear behavior, we can write the piezoelectric constitutive equation as:


The Dirichlet and the Neumann boundary conditions for the mechanical field, see Fig. 4, are


where and are prescribed mechanical displacement and surface traction vectors in the reference configuration, and , , with as the boundary of the domain and , as its Dirichlet and Neumann portions. Moreover N is the outward unit normal to . The boundary conditions for the electric field are


where and are prescribed values of electric potential and electric charge flux, and , .

3.2 Finite element formulation including electromechanical contact

Within the RVE, the piezoelectric fibers are subjected to contact constraints. For the sake of simplicity, we assume herein frictionless contact conditions. Since no experimental results are available on contact interactions between the piezoelectric fibers examined in this paper, the choice of a friction law and of the related coefficient(s) would have been arbitrary. The addition of frictional effects and the corresponding experimental justification may well be considered in further developments of this research. Herein the main characteristics of the 3D electromechanical frictionless contact formulation are now summarized. All contact related quantities are referred to the micro scale. The formulation is based on the classical master-slave concept [30], see Fig. 5. For each point on the slave surface, , the corresponding point on the master surface is determined through normal (i.e. closest point) projection and is denoted as . Thus the normal gap of each slave point is computed as:


being the outer unit normal to the master surface at the projection point in the current configuration. 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 and are the electric potential values in the slave point and in its master projection point.

The discretization strategy used herein for the contact contribution is based on the node-to-surface approach combined with B zier smoothing of the master surface. It is well known that the node-to-surface algorithm is susceptible of pathologies due to the -continuity of the finite element (Lagrange) discretizations, and that these may affect the quality of results as well as iterative convergence in contact computations [7, 32]. One of the possible remedies are smoothing techniques for the master surface. Herein, the technique based on B zier patches proposed by [17] and implemented within the AceGen/AceFEM environment is adopted and straightforwardly extended to electromechanical contact constraints, see Fig. 5. Here a tensor product representation of one-dimensional B zier polynomials is used to interpolate the master surfaces in the contact interface. A three-dimensional B zier surface of order is given by


where and are Bernstein polynomials, functions of the convective surface coordinates , and are the coordinates of the so-called control points. To fully characterize a bicubic B zier patch, 16 control points are needed.

The technique used herein consists in the so-called B zier-9 patches [17]. Each patch is constructed using one central node of the master surface, denoted as , and eight neighbouring nodes , () forming four bilinear segments surrounding the central node. The 16 control points needed to define the bicubic B zier patch are obtained from


where , is a parameter defining the shape of the surface (here ), and the auxiliary points are defined in terms of the master nodes according to


The same procedure is followed for the electric potential, which is expressed as


where now is the potential evaluated at the 16 control points as a function of the potential at the auxiliary points and at the master nodes .

Note that the need for smoothing of the master surface can be straightforwardly eliminated by the use of isogeometric discretizations, as these feature higher continuity and smoothness than classical Lagrange discretizations (see e.g. [8]). This will be done in future extensions of the present work.

According to standard finite element techniques, the global set of equations can be obtained by adding to the variation of the energy potential representing the continuum behavior the virtual work due to the electromechanical contact contribution associated to the active contact elements. The global energy of the discretized system (where refers again to the microscale) is thus


The virtual work due to contact is in turn given by


with and as the virtual work due to mechanical and electric contact, respectively. Herein, the electromechanical contact constraints are regularized with the penalty method, which leads to




are respectively the contact force and the electric current, and and are penalty parameters. Virtual variation of eq. (49) gives


and its consistent linearization leads to


The advanced symbolic computational tools available in the AceGen/AceFEM finite element environment allow for a full automation of the linearization process, see [14] for more details. If and are the discretized displacement and electric potential fields, with and as the nodal displacements and electric potential, respectively, and as the (in this case linear) shape functions, the residual vector and the stiffness matrix terms resulting from the finite element discretization are determined as follows


4 Computational homogenization procedure

In this section, a two-step homogenization method is introduced. Starting from the microscale RVE, we obtain in the first step the macroscopic constitutive response of an equivalent homogenized solid, and in the second step, through thickness integration, the effective coefficients of an homogenized shell. First, the concept of RVE and Hill’s energy principle are briefly introduced. This allows for the transition from micro to macro quantities characterizing the multiphysics problem. Then, a method to calculate the macroscopic constitutive matrix starting from the solution of a microscale boundary value problem (BVP) is described for the general case of electromechanical problems.

4.1 Unit cell models for numerical homogenization

The main idea of homogenization is finding an homogeneous medium equivalent to the original heterogeneous material, such that the strain energies stored in the two systems are the same (Hill condition). Split of the microscopic deformation gradient and electric field into constant parts and fluctuating parts yields:


where by definition


In the previous and in the next equations (unless otherwise specified), the RVE volume in the reference configuration , Fig. 6, is indicated as for the sake of a compact notation. Formulated for the electromechanical problem at hand, the Hill 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 the work on the macroscale. In the previous equation as well as in the following, the superscript M refers to the macroscale. In particular, , , , and represent respectively the average values of the first Piola Kirchhoff stress, electric displacement, deformation gradient and electric field, i.e.


Eq. 64 may be split into two parts


4.2 Boundary conditions

Classically three types of boundary conditions are used for an RVE [2, 26]: prescribed linear displacements, prescribed constant tractions and periodic boundary conditions. All three types of boundary conditions satisfy the micro-macro work equality stemming from Hill’s lemma in eq. (64) and are therefore suitable for the analysis. Periodic boundary conditions (Fig. 7) are suitable for solids with periodic microstructure (Fig. 2), such as the PVDF sheet under study, and therefore are preferred in this paper. The periodicity conditions for the RVE are derived in a general format as:


where the indices and indicate the values on two opposite surfaces of the unit cell, and , , , , , , and , are respectively the components of displacement, traction, electric potential, and electric flux at corresponding points on the opposing faces and (Fig. 7). Periodic boundary conditions are implemented in Acegen as multipoint constraints using the Lagrange multiplier method. Uniform meshing of the RVE facilitates the enforcement of these constraints, as opposing nodes match on all the faces of the RVE. Note that the RVE structure in this study may be considered as a particular case of bimaterial RVE, where the first material corresponds to the fibers () and the second material is represented by the void between the fibers in contact (). In Fig. 6, the total RVE volume is thus . During volume integration the void part affects the average values of the electromechanical properties. Perfect bond between the fibers is assumed along the generatrices of the cylinders (blue lines in Fig. 6) and is enforced through condensation of the respective degrees of freedom. This assumption is based on observations of the fibers with the Scanning Electron Microscope, which reveal as “welding” takes place between the fibers upon manufacturing due to solvent evaporation [22].

4.3 Macroscale constitutive law

Due to the electromechanical contact between fibers at the microscale, the macroscopic constitutive law is nonlinear but can be linearized introducing a tangent or secant constitutive matrix for the macroproblem. In this work, we opted for the secant approach, see also [37, 38, 39], in which the effective (or volume-averaged) generalized strain column vector of the homogenized material is computed successively (and independently) for increasing values of the effective generalized stress vector . They are related in each step through a generalized secant effective constitutive matrix so that




where , , and are the matrices of the secant elastic, piezoelectric and dielectric coefficients. The internal energy function per unit volume for the solid is:


where the superscript stands for solid. If now indicates the total shell volume , the potential of the internal forces acting on the shell can be obtained through integration as:


where eq. (25) has been used and, in the last equality, the macroscopic constitutive matrix of the shell has been defined as


where is the determinant of the shifter tensor that for slightly curved shell is . On the other hand, the potential of the shell can be expressed as


where is the vector of stress resultants of the shell given by


where are the membrane forces, the bending moments, the shear forces, the electric displacements, whereas , are the constant and linear parts of the components in the thickness direction. Therefore, the comparison between eqs. (73) and (75) leads to the expression of as


5 Results

As follows, the presented computational procedure is applied to a simple test case, whereby an RVE containing four PVDF fibers in a square configuration is considered. The final aim is to predict the effective properties of the homogenized shell element. Each fiber is considered as a linear piezoelastic solid and is discretized with linear 8-node brick elements (Fig. 9). Similarly to [9] the nonlinear theory is specialized to the kinematically linear case where the mechanical constitutive behavior is defined by the following constitutive law:


where denotes the Cauchy stress, the linearized strain tensor with trace . Moreover E is the Young’s modulus and is the Poisson’s ratio. Frictionless electromechanical contact constraints are enforced at the interface between the fibers using the described electromechanical contact formulation. The final piezoelectric constitutive matrix of the PVDF material, , simplifies to