# Electromechanical homogenization of polyvinylidene fluoride PVDF polymer reinforced by graphene nanoribbons

Salah Elbarnaty, Wiyao Azoti, Joao Pedro M Correia, Nadia Bahlouli, Saïd Ahzi

► To cite this version:

Salah Elbarnaty, Wiyao Azoti, Joao Pedro M Correia, Nadia Bahlouli, Saïd Ahzi. Electromechanical homogenization of polyvinylidene fluoride PVDF polymer reinforced by graphene nanoribbons. Modelling and Simulation in Materials Science and Engineering, 2025, 33 (6), pp.065011. 10.1088/1361-651X/adf870. hal-05280924

HAL Id: hal-05280924

https://hal.science/hal-05280924v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**PAPER • OPEN ACCESS**

Electromechanical homogenization of
polyvinylidene fluoride PVDF polymer reinforced
by graphene nanoribbons

To cite this article: Salah Elbarnaty et al 2025 *Modelling Simul. Mater. Sci. Eng.* **33** 065011

View the [article online](#) for updates and enhancements.

You may also like

- Photonic-digital hybrid artificial intelligence hardware architectures: at the interface of the real and virtual worlds
  Lilia M S Dias, Dinis O Abranches, Ana R Bastos et al.

- ICRH modelling of DTT in full power and reduced-field plasma scenarios using full wave codes
  A Cardinali, C Castaldo, F Napoli et al.

- Global evidence that cold rocky landforms support icy springs in warming mountains
  Stefano Brighenti, Constance I Millar, Scott Hotaling et al.
# Electromechanical homogenization of polyvinylidene fluoride PVDF polymer reinforced by graphene nanoribbons

Salah Elbarnaty<fnref n="1" />,* , Wiyao Azoti<fnref n="1" />,<fnref n="2" />,
Joao Pedro M Correia<fnref n="1" /> , Nadia Bahlouli<fnref n="1" /> and Saïd Ahzi<fnref n="1" />,<fnref n="3" />

E-mail: s.elbarnaty@unistra.fr

Received 26 May 2025; revised 23 July 2025

Accepted for publication 6 August 2025

Published 14 August 2025

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Our work consists in modeling the effective properties of a piezoelectric-based nanocomposite that will be used for aeronautical applications. For this purpose, the polyvinylidene fluoride (PVDF) is considered as a polymer piezoelectric matrix in which graphene nanoribbons (GnR) are spatially oriented as inclu- sions. Two different methods were used to estimate the effective properties: the Mori-Tanaka (MT) micromechanics scheme and the finite element (FE) ana- lysis. In the MT homogenization scheme, the electro-elastic Eshelby tensor was first computed for an anisotropic matrix. The entire MT scheme was implemen- ted in the MATLAB environment. Afterwards a FE model of a nanocomposite unit cell was developed with ABAQUS in order to estimate the effective prop- erties. The effective properties computed with the MT scheme showed good agreement with the predictions of the FE model. The orientation and volume fraction of GnR inclusions in the PVDF matrix were found to play a signific- ant role in enhancing the mechanical properties of the nanocomposite, with a moderate effect on its piezoelectric response. In addition, several physical mechanisms such as stacking-up and alignment defects of GnR were analyzed

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ ICUBE laboratory-CNRS UMR 7357, University of Strasbourg, CS 90032, Strasbourg 67081, France</footnote>

<footnote n="2">² Now at INSA Toulouse, ICA, CNRS UMR 5312, Toulouse, France</footnote>

<footnote n="3">³ LERMA laboratory, International University of Rabat, Sala Al Jadida 11100, Morocco</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
for their negative influence on the reinforcement effect of GnR. Our theoretical study will help in understanding the favorable effects of nano-inclusions over piezoelectric nanocomposites.

Keywords: piezoelectricity, homogenization, FE analysis, Mori–Tanaka method, PVDF, GnR

# 1. Introduction

For the last decades interest in smart materials is on the increase. Among them, piezoelectric materials are becoming sought after because they are used for producing various devices such as sensors or actuators in applications such as medical imaging scanners, electronic devices and wearable technology, their scope of application has recently expanded significantly. For instance, piezoelectric materials have attracted increasing interest in recent years in aeronautical and space applications [1]. In works by Li et al [2] it was theoretically demonstrated that piezoelectric actuators could reduce the aircraft's roll motion while saving structure weight. Furthermore, piezoelectric materials are essential in the development of the morphing wing concept [3, 4]. More recently, numerical works by Abou-Khalil et al [5] have demonstrated that the aerodynamic performance of an aircraft could be improved by a moving interface between the aircraft wing profile and the surrounding turbulence. The moving or vibrating interface would consist in a large number of electro-active tiles made of a new piezoelectric material with high sensing and actuation capacities. To limit the weight of the interface, it should be as thin as possible and made of a piezoelectric nanocomposite with a polymer matrix reinforced by nanoinclusions (in order to improve the piezoelectric properties).

In this context, polyvinylidene fluoride (PVDF) emerges as a notable alternative to piezoceramic material (e.g. PZT family), mainly because the latter is lead-based material [6]. Beyond its good piezoelectric properties [7, 8], PVDF is thermally stable and has favorable mechanical characteristics [8]. In its natural state, PVDF exists in an amorphous phase (denoted α phase) [9], requiring a process to transform it into an electro-active state (crystal-line phase such as β phase)—the phase of particular interest for piezoelectricity [10]. Various methods, such as mechanical stretching, high voltage polarization or thermochemical treatments can be used to activate the transformation phase process [9–11]. However, PVDF can also be obtained by electrospinning to enhance β phase formation [8].

The mechanical and piezoelectric properties of PVDF can be significantly improved by adding nanoinclusions such as carbon nanotubes (CNT) or graphene nanoplatelets (GnP). Several experimental or numerical studies have highlighted the effects of nanoinclusions on the behaviour of PVDF [7, 12–14]. Besides the significant and obvious mechanical improvement, the piezoelectric coefficient $d_{33}$ could be increased from 27.4 pC N⁻¹ to 31.3 pC N⁻¹ by just adding 0.01 Wt% of CNT [15]. Nevertheless, nanoinclusions have high electrical conductivity [16] and therefore must be added in very small amounts due to the percolation obstacle [17, 18], which could transform the piezoelectric polymer matrix into a conductive material. More recently, researchers have attempted to overcome the problem by creating defects in the nanoinclusions in order to transform them into dielectric or even piezoelectric materials [19–21]. For example, graphene nanoribbons (GnR), obtained by unzipping CNT, have high dielectric properties [21, 22] and thus could be used as nanoinclusions for our intended application: a vibrating interface made of a piezoelectric nanocomposite.

In order to design a new piezoelectric nanocomposite, it is necessary to establish accurate structure-property expressions. Extensive efforts have been made to estimate the effective properties of piezoelectric composites. The first approaches use homogenization methods
based on classical Eshelby's inclusion solution and numerical schemes to compute the effect-
ive properties [23, 24]. Analytical relationships have also been developed such as those by
Kuo and Huang [25]; they show a strong coupling between the orientation and the volume
fraction of inclusions with the effective piezoelectric properties. Expressions for Eshelby's
electro-elastic tensor are proposed in [26, 27] and used to derive the self-consistent model, the
Mori-Tanaka (MT) scheme and the dilute approach.

According to previous works [23, 26, 28], the MT averaging micromechanics scheme,
which was initially developed to deal with the mechanical behavior of composites and was
later extended to include electro-mechanical coupling, shows good agreement with experi-
mental data. Odegard et al [29] employed the MT method to predict the mechanical behavior
of effective CNT fibers, which were first characterized using molecular dynamics simulations
for both aligned and randomly distributed inclusions. Their predictions showed good agree-
ment with experimental data, particularly at low volume fractions of randomly distributed
CNTs. Ji et al [30] applied the MT scheme to predict the effective properties of GnP-based
composites with randomly distributed inclusions. Similarly, Azrar et al [31] and Li and Dunn
[32]. utilized the MT method to model the viscoelectroelastic behavior of two-phase compos-
ites. Furthermore, Li et al [33] employed the MT approach to predict the viscoelastic behavior
of CNT-reinforced polymer composites in the frequency domain (Carson domain).

By employing the approach of representative volume element (RVE), numerical methods
such as finite element (FE) or boundary element methods have been used to estimate the effect-
ive properties of nano composites and piezoelectric composites [28, 34–37]. According to
the literature, the FE method is also found to be accurate because it takes into account the
stress distribution within the inclusion, unlike Eshelby-based schemes which assume stress
uniformity in the inclusion. Chwał and Muc [35, 38] applied FE micromechanical model-
ing to characterize the elastic behavior of CNT polymer nanocomposites, emphasizing the
impact of nanotube arrangement. Kumar and Srivastava [39] used a 3D nanoscale RVE with
FE approach to evaluate elastic properties of CNT- and GnR-reinforced composites, analyzing
effects of volume fraction and reinforcement length and compared them with continuum mech-
anics models (i.e. Halpin–Tsai and the rule of mixtures models). Guo et al [36] carried out a
multiscale FE model combining atomistic and macroscopic scales, and effectively predicted
that increasing graphene volume fraction and size significantly enhance the Young’s modulus
and stress transfer in graphene-reinforced composites. El Bahi et al [40] introduced a multi-
stage homogenization strategy to assess how CNT/GnP agglomeration and geometry influence
the elastic properties of polymer composites reinforced with both nano- and micro-scale fillers.
Shingar and Naskar [41] investigated the FE approach and various micromechanical modeling
techniques (i.e. Rule of mixtures, modified rule of mixtures and Halpin–Tsai approaches) to
predict the enhanced elastic and piezoelectric properties of graphene-reinforced epoxy com-
posites and PZT-based composites, demonstrating that graphene nanofillers—through their
shape, orientation, and strong interfacial bonding—significantly improve composite perform-
ance, with results closely matching experimental data. Most recently, FE approach has gained
increased attention as a reliable reference method for generating training data in machine learn-
ing models designed to predict the mechanical and piezoelectric behavior of composites [42–
46]. Together, these studies validate the Finite Element method as a powerful and accurate tool
to assess nanocomposite properties.

Piezoelectric GnR/PVDF nanocomposite materials are studied in our work. First, a brief
introduction on the constitutive law of linear piezoelectricity is developed. Then, a homogen-
ization modeling based on the MT micromechanics scheme is presented in order to predict the
effective properties of piezoelectric nanocomposites. For comparison purposes, 3D FE models
are also developed using the commercial FE code ABAQUS/CAE 2020. The effects of orientation and volume fraction of GnR inclusions are studied as well as GnR stacking and alignment defects in terms of change in the effective piezoelectric properties. The present study is limited to nanocomposites without defects and with perfectly straight GnR inclusions. Moreover, all GnR inclusions are assumed to be perfectly embedded into the PVDF matrix.

## 2. Materials

The studied piezoelectric nanocomposites consist of a PVDF matrix reinforced by nano-inclusions of graphene. Graphene inclusions can have different geometries such as nanoplatelets (2D geometry and called GnP) or nanoribbons (quasi-1D geometry called GnR). In the case of nanoribbons, graphene appears in the form of narrow and very long strips. The mechanical properties of GnR are generally evaluated from numerical simulations (using molecular dynamics or finite element analysis). GnR behave as an anisotropic elastic material with elastic properties that depend on its geometry (length and atomic structure) [47]. However, the Young's modulus and Poisson ratio seem to converge as the length of GnR increases [47]. For simplicity, GnR are assumed to have an isotropic linear elastic behavior with constant properties. PVDF is a piezoelectric polymer with a semi-crystalline structure and is assumed to have a large fraction of electro-active β-phase. Moreover, PVDF shows a nonlinear behavior that is time and temperature dependent [48–50]. However, at small deformations, its mechanical behavior can be considered linear elastic.

In the work, PVDF is considered as a continuous matrix in which GnR inclusions are embedded with a small volume fraction. The latter assumption is justified because the GnR volume fraction must not reach the agglomeration and percolation threshold constraint [51]. It is also assumed that the addition of GnR inclusions does not change the fraction of electro-active β-phase of PVDF. Therefore, the matrix material properties are kept constant. The properties of both PVDF and GnR are taken from the literature [28, 40, 52], and are listed in table 1.

## 3. Constitutive equations for piezoelectric materials

The constitutive equations of piezoelectric materials have a coupled interaction between the electrical and mechanical behaviors. For a linear piezoelectric material, the local constitutive behavior is given in the tensorial form by [26, 28]:

$$ 
\begin{aligned}
\sigma_{ij} &= C_{ijmn}\varepsilon_{mn} + e_{nij}E_n \\
D_i &= e_{imn}\varepsilon_{mn} - \kappa_{in}E_n
\end{aligned}
\qquad (1) $$

where $\sigma_{ij}$, $\varepsilon_{mn}$, $E_n$, and $D_i$ are the 2nd order stress tensor, the 2nd order strain tensor, the 1st order electric tensor and 1st order electric displacement tensor respectively. Whereas $C_{ijmn}$, $e_{ijk}$ and $\kappa_{in}$ are 4th order elastic tensor, the 3rd order piezoelectric tensor and 2nd order dielectric tensor, respectively.

In the absence of other external forces, the divergence of local constitutive behavior leads to the equilibrium equations:

$$ 
\begin{aligned}
\sigma_{ij,j} &= 0 \\
D_{i,i} &= 0.
\end{aligned}
\qquad (2) $$
<table><caption>**Table 1.** Elastic, piezoelectric, and dielectric properties of the PVDF and GnR ($k_0 = 8.85$ F m$^{-1}$, electric permittivity of free space).</caption><thead><tr><th>Properties</th><th>PVDF [28]</th><th>GnR [40, 53]</th></tr></thead><tbody><tr><td>C<sub>11</sub> (GPa)</td><td>3.8</td><td>1600.0</td></tr><tr><td>C<sub>12</sub> (GPa)</td><td>1.9</td><td>864.2</td></tr><tr><td>C<sub>13</sub> (GPa)</td><td>1.0</td><td>864.2</td></tr><tr><td>C<sub>22</sub> (GPa)</td><td>3.2</td><td>1600.0</td></tr><tr><td>C<sub>23</sub> (GPa)</td><td>0.9</td><td>864.2</td></tr><tr><td>C<sub>33</sub> (GPa)</td><td>1.2</td><td>1600.0</td></tr><tr><td>C<sub>44</sub> (GPa)</td><td>0.7</td><td>370.0</td></tr><tr><td>C<sub>55</sub> (GPa)</td><td>0.9</td><td>370.0</td></tr><tr><td>C<sub>66</sub> (GPa)</td><td>0.9</td><td>370.0</td></tr><tr><td>k<sub>1</sub>/k<sub>0</sub></td><td>7.4</td><td>6.9</td></tr><tr><td>k<sub>2</sub>/k<sub>0</sub></td><td>9.3</td><td>6.9</td></tr><tr><td>k<sub>3</sub>/k<sub>0</sub></td><td>7.6</td><td>6.9</td></tr><tr><td>e<sub>31</sub> (C m<sup>-2</sup>)</td><td>0.024</td><td>0.0</td></tr><tr><td>e<sub>32</sub> (C m<sup>-2</sup>)</td><td>0.001</td><td>0.0</td></tr><tr><td>e<sub>33</sub> (C m<sup>-2</sup>)</td><td>-0.027</td><td>0.0</td></tr></tbody></table>

The equations for the elastic strain and the electric field are given, respectively, as follows:

$$ \varepsilon_{ij} = \frac{1}{2} (u_{i,j} + u_{j,i}) \qquad (3) $$

$$ E_i = -\psi_{,i} $$

where $u_i$ and $\psi$ are the mechanical displacement and the electric potential. Following conventional subscript notation, lower-case subscripts are considered from 1–3, while upper-case subscripts take range from 1–4, and repeated lower-case (or upper-case) subscripts sum over 1–3 (or 1–4). Under this convention, the electric and elastic variables are combined to form a compact linear equation such as:

$$ \Sigma_{iJ} = E_{iJMn}Z_{Mn} \qquad (4) $$

where $Z_{Mn}$ states for the elastic-strain and electric fields

$$ Z_{Mn} = \begin{cases} \varepsilon_{mn} & M = 1, 2, 3 \\ -E_n & M = 4 \end{cases} \qquad (5) $$

while $\Sigma_{iJ}$ represents the elastic stress and electric displacement such as:

$$ \Sigma_{iJ} = \begin{cases} \sigma_{ij} & J = 1, 2, 3 \\ D_i & J = 4 \end{cases} \qquad (6) $$

The electro-elastic moduli $E_{iJMn}$ are given by:

$$ E_{iJMn} = \begin{cases} C_{ijmn} & J = 1, 2, 3; \quad M = 1, 2, 3 \\ e_{nij} & J = 1, 2, 3; \quad M = 4 \\ e_{imn} & J = 4; \quad M = 1, 2, 3 \\ -\kappa_{in} & J = 4; \quad M = 4 \end{cases} \qquad (7) $$

The equation (4) can be reduced to a matrix form using the Voigt notation which gives:

$$ \Sigma = EZ. \qquad (8) $$
The bold character indicates a matrix form and consequently, **E** stands for a (9 × 9) matrix while **Σ** and **Z** are (9 × 1) vectors. The final matrix system for a linear orthotropic piezoelectric material can be written as:

$$ \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \\ D_1 \\ D_2 \\ D_3 \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 & 0 & 0 & e_{31} \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 & 0 & 0 & e_{32} \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 & 0 & 0 & e_{33} \\ 0 & 0 & 0 & C_{44} & 0 & 0 & 0 & e_{15} & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 & e_{15} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & e_{15} & 0 & 0 & 0 & 0 \\ e_{31} & e_{32} & e_{33} & 0 & 0 & 0 & 0 & 0 & -\kappa_3 \end{bmatrix} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \\ E_1 \\ E_2 \\ E_3 \end{bmatrix} \quad (9) $$

where $2\varepsilon_{ij} = \gamma_{ij}$, $i \neq j$,

# 4. Effective electro-elastic properties

## 4.1. MT micromechanics scheme

The effective properties are derived from the local constitutive equations along with a strain localization step through the strain concentration tensor. Using an average technique over the local equations within a RVE, the effective properties based on the MT scheme are obtained explicitly using the equivalent inclusion of Eshelby [54]. For electro-elastic composites, details about the development of the MT scheme by volume averaging can be found in the literature [23, 26]. In our study, the effective properties using the MT scheme for a two-phase composite are given as:

$$ \mathbf{E}^{\mathrm{MT}} = \mathbf{E}_{\mathrm{PVDF}} + \mathbf{c}_{\mathrm{GnR}} (\mathbf{E}_{\mathrm{GnR}} - \mathbf{E}_{\mathrm{PVDF}}) \mathbf{A}^{\mathrm{MT}} \quad (10) $$

with

$$ \mathbf{A}^{\mathrm{MT}} = \mathbf{A}^{\mathrm{dil}} \left[ \mathbf{c}_{\mathrm{PVDF}} \mathbf{I} + \mathbf{c}_{\mathrm{GnR}} \mathbf{A}^{\mathrm{dil}} \right]^{-1} \quad (11) $$

where $\mathbf{c}_{\mathrm{PVDF}}$, $\mathbf{c}_{\mathrm{GnR}}$ are the matrix and inclusion volume fractions respectively. The matrices $\mathbf{E}_{\mathrm{PVDF}}$ and $\mathbf{E}_{\mathrm{GnR}}$ are respectively the electro-elastic matrices of the PVDF and the GnR. $\mathbf{E}^{\mathrm{MT}}$ is the effective electromechanical matrix predicted by the MT method. In equation (10), $\mathbf{A}^{\mathrm{MT}}$ is the strain concentration matrix. The matrix $\mathbf{A}^{\mathrm{dil}}$ is the strain concentration matrix when the inclusions are diluted in the matrix without interaction towards each other [55] and it is given by:

$$ \mathbf{A}^{\mathrm{dil}} = \left[ \mathbf{I} + \mathbf{SE}_{\mathrm{PVDF}}^{-1} (\mathbf{E}_{\mathrm{GnR}} - \mathbf{E}_{\mathrm{PVDF}}) \right]^{-1} \quad (12) $$

where **I** is the (9 × 9) identity matrix while **S** is the matrix form of the 4th-order piezoelectric Eshelby tensor. Its tensorial expression $S_{MnAb}$ is given by the equation (13). It depends on the matrix electro-elastic properties and on the geometry of the inclusion. The $S_{MnAb}$ components for an ellipsoidal inclusion are given in [23, 24] such as:

$$ S_{MnAb} = \begin{cases} \frac{a_1 a_2 a_3}{8\pi} E_{iJAb} \int_{|z|=1} \frac{1}{\zeta^3} [G_{mJin}(z) + G_{nJin}(z)] dS(z), & M = 1, 2, 3 \\ \frac{a_1 a_2 a_3}{4\pi} E_{iJAb} \int_{|z|=1} \frac{1}{\zeta^3} G_{4Jin}(z) dS(z), & M = 4 \end{cases} \quad (13) $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Typical scheme for an ellipsoidal inclusion embedded in an infinite matrix.</figcaption></figure>

where $a_1, a_2$ and $a_3$ are the half axes of the ellipsoidal inclusion and $|z| = 1$ is the surface of the unit sphere.

$$ \zeta = [a_1^2 z_1^2 + a_2^2 z_2^2 + a_3^2 z_3^2]^{1/2} \quad (14) $$

where $G$ stands for the Green tensor, which is given by:

$$ G_{MSin}(z) = z_i z_n K_{MJ}^{-1}(z) \quad (15) $$

with:

$$ K_{MJ}(z) = E_{iJMn} z_i z_n. \quad (16) $$

The unit sphere is parameterized by $\xi_3$ and $\theta$ as follows:

$$ z_1 = \frac{\xi_1}{a_1} = \frac{(1 - \xi_3^2)^{1/2} \cos\theta}{a_1}, z_2 = \frac{\xi_2}{a_2} = \frac{(1 - \xi_3^2)^{1/2} \sin\theta}{a_2}, z_3 = \frac{\xi_3}{a_3}. \quad (17) $$

In our work the components of $S$ tensors are computed numerically since the PVDF is an orthotropic material. It should be noted that an explicit analytical solution exists only in the case of a transversely isotropic matrix [5]. In general, the tensor $S$ is not symmetric and care must be taken when it comes to the shear part [23].

In our study, the GnR are considered as nanoribbons with $a_2 \gg a_1$ and $a_3 = \infty$. GnR are spatially oriented within the PVDF Matrix. Two coordinate systems are considered. The first one is associated to the inclusion and is noted ($x'_1, x'_2, x'_3$). The second one is the global coordinate system. The latter is related to the whole composite and is called ($x_1, x_2, x_3$) as shown in figure 1. The non zero-elements of the Eshelby matrix of GnR reinforced PVDF composite can be found in appendix A.
The two coordinate systems are connected by the Euler angles $\phi$, $\phi$ and $\omega$ with a 323-transformation type such as:

$$ \begin{pmatrix} x'_{1} \\ x'_{2} \\ x'_{3} \end{pmatrix} = [\alpha_{ij}] \begin{pmatrix} x_{1} \\ x_{2} \\ x_{3} \end{pmatrix} \qquad (18) $$

where $[\alpha_{ij}]$ is the cosine direction matrix of type 323 and is given by:

$$ [\alpha_{ij}] = \begin{bmatrix} \cos \varphi \cos \phi \cos \omega - \sin \varphi \sin \varphi & \sin \varphi \cos \phi \cos \omega + \cos \varphi \sin \varphi & -\sin \phi \cos \omega \\ -\cos \varphi \cos \phi \sin \varphi - \sin \varphi \cos \omega & \cos \varphi \cos \omega - \sin \varphi \cos \phi \sin \varphi & \sin \phi \sin \varphi \\ \cos \varphi \sin \phi & \sin \varphi \sin \phi & \cos \phi \end{bmatrix}. \quad (19) $$

It is worth noting that the transformation described above should not be applied directly on $E_{iJMn}$ nor $S_{MnAb}$, because none of them is a tensor in the true sense. In the case of $E_{iJMn}$, the transformation is applied directly on each tensor of the set defined by equation (20) (because of the lack of the major symmetry). Subsequently, the Voigt notation is used to convert each tensor into the matrix form. The transformation is applied on the elastic, piezoelectric, and dielectric parts of the electro-elastic (given an inclusion orientation) as follows:

$$ C'_{ijkl} = C_{mnpq} \alpha_{im} \alpha_{jn} \alpha_{kp} \alpha_{lq} \qquad (20) $$

$$ e'_{ijk} = e_{mnp} \alpha_{im} \alpha_{jn} \alpha_{kp} $$

$$ \kappa'_{ij} = \kappa_{mn} \alpha_{im} \alpha_{jn}. $$

Finally, all the calculations discussed in section 4 were processed with MATLAB using the algorithm given by figure 2.

## 4.2. FE modeling

In FE simulations, the strain and stress fields are computed in the inclusion and also in the matrix, which leads to more accurate predictions of the effective properties [28]. Therefore, FE predictions are used to validate the results of the micro-mechanical modeling method presented in the previous section. The 3D FE model is established using the commercial FE code ABAQUS/CAE 2020.

The RVE of the nanocomposite is built in accordance with [39] considering a perfectly straight GnR. Its 3D geometric configuration is illustrated in figure 3. In the RVE, the GnR is perfectly embedded into the PVDF matrix (perfect interface between them) and arranged along its poling direction (3-direction). The interphase zone between the inclusion and the matrix is neglected in FE modeling. The change in volume fraction is made by modifying the size of the surrounding matrix and keeping the size of the GnR constant. The RVE is meshed with 8-node linear piezoelectric brick elements (called C3D8E in ABAQUS/CAE element library). An example of the RVE mesh is presented in figure 4. The RVE mesh corresponds to a GNR volume fraction equal to 2.45% which is discretized into 612 200 elements. The materials properties inputs used in the FE analysis are listed in table 1.

Various loading conditions have to be considered to estimate the effective properties of the nanocomposite. The boundary conditions corresponding to the loading conditions must be carefully defined and are derived from [28]. Two types of boundary conditions are imposed on the symmetric faces of the FE model: mechanical displacement $u_i$ and electric potential $\psi$. They are respectively defined by:

$$ u_i(F) = \varepsilon_{ij}^0 x_j \qquad (21) $$

$$ \psi(F) = -E_j^0 x_j $$
<figure><img src="image_5.png" /><figcaption>Figure 2. Algorithm for computing the MT effective electro-elastic properties.</figcaption></figure>

where corresponds to the face of the FE mesh on which the boundary conditions are imposed, $ε_{ij}^0$ the constant applied strain, $E_j^0$ the constant applied electric field and $x_j$ the coordinates of the nodes belonging to the RVE's face F.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Geometric configuration of the RVE showing an inclusion of GnR embedded in the PVDF matrix.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** FE mesh of the RVE in the case of a GnR/PVDF nanocomposite with a GnR volume fraction equal to 2.45%.</figcaption></figure>

The effective properties of the nanocomposite are calculated by equating the numerical predictions of total energies with their analytical expressions. In the case of effective elastic properties, the elastic energy is first determined for each element of the FE mesh of the RVE. Then, by summing over all the RVE elements, the total elastic energy stored in the RVE $U_e$ is deduced as shown in equation (22). Finally, the equivalent elastic properties are obtained by solving the equality with the analytical expression, equation (22). The effective piezoelectric
properties are similarly derived from the total electro-mechanical energy stored in the RVE $U_{em}$, equation (23). The total elastic and electro-mechanical energies are defined by:

$$U_e = \sum_{m=1}^{n} U_e^m = (V/2) C_{ijkl} \varepsilon_{ij}^0 \varepsilon_{kl}^0 \quad (22)$$

$$U_{em} = \sum_{m=1}^{n} U_{em}^m = (V/2) e_{ijk} \varepsilon_{jk}^0 E_i^0 \quad (23)$$

where the index $m$ corresponds to the finite element considered in the FE mesh, the total number of finite elements and $V$ the total volume of the RVE. In ABAQUS/CAE 2020, $U_{em}^m$ is the sum of the elastic strain energy and the dielectric energy stored in the element $m$ of the RVE. In ABAQUS/CAE documentation, they are referred to as ELSE and EENER respectively. A typical Python code used for postprocessing in order to obtain the effective electro-elastic constants is given in the appendix B. Tables 2 and 3 summarize the boundary conditions and associated formulas used to estimate the effective elastic and piezoelectric constants by FE analysis.

# 5. Results and discussion

## 5.1. Comparison between MT and FE methods

There is a lack of experimental results regarding the effective properties of GnR/PVDF nano-composites in the literature. In order to check the validity of the developed micromechanical approach, the electro-elastic properties of the nanocomposite GnR/PVDF evaluated by the MT method are compared with those evaluated by the FE method. Consequently, the effective elastic moduli $Y_{11}$, $Y_{22}$, $Y_{33}$, $G_{12}$, $G_{13}$ and $G_{23}$, as well as the effective piezoelectric constants $e_{31}$, $e_{32}$ and $e_{33}$, are presented as a function of the volume fraction of GnR in figures 5 and 6. The configuration in figures 5 and 6 corresponds to all GnR inclusions perfectly aligned with the third material direction of the PVDF matrix (i.e. the poling direction of PVDF), a configuration was chosen in order to enhance the piezoelectric effect in lateral direction (the direction of interest for actuation mode). In the MT method, the alignment of GnR inclusions implies that the Euler angles ($\varphi$, $\phi$ and $\omega$) are set to 0 and therefore that the matrix $[a_{ij}]$ is reduced to the identity matrix in equation (21).

In figures 5 and 6, the predictions of all electro-elastic properties computed by the MT method show good agreement with those evaluated by FE simulations. The proposed micromechanical modeling based on the MT scheme predicts an almost constant trend of the piezoelectric constants $e_{ij}$ as the volume fraction of GnR increases, figure 6, whereas it predicts a strong sensitivity of Young's and shear moduli. As expected, the effective values of Young's moduli linearly increase when the volume fraction of GnR inclusions increases, figure 5(a). The effective values of shear moduli also linearly increase with increasing GnR volume fraction, figure 5(b). The linear increase of the effective elastic properties due to the increase in the graphene volume fraction is in agreement with some numerical works in the literature, such as in [7]. The strong improvement of Young's and shear moduli is due to the high mechanical properties contrast of PVDF and GnR. From figures 5(a) and (b), it can also be concluded that the elastic properties of the nanocomposite become highly anisotropic. The effective piezoelectric constants $e_{ij}$ predicted by both methods have fairly constant values, figure 6. They are therefore not altered by the addition of GnR. These numerical predictions can be explained: in the proposed modeling it was assumed that the $\beta$-phase fraction of the PVDF matrix was not modified by the GnR volume fraction.
<table><caption>**Table 2.** Boundary conditions for elastic constants [28]. The +and−respectively refer to the front face and back face of the RVE for a given specific direction. The electric fields are taken equal to zero.</caption><thead><tr><th>Elastic Constants</th><th>Applied deformation</th><th>Applied displacement</th><th>Elastic energies</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>ε<sup>0</sup><sub>11</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = ε<sup>0</sup>x<sub>1</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>11</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>22</sub></td><td>ε<sup>0</sup><sub>22</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = ε<sup>0</sup>x<sub>2</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>22</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>33</sub></td><td>ε<sup>0</sup><sub>33</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = ε<sup>0</sup>x<sub>3</sub>, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>33</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>12</sub></td><td>ε<sup>0</sup><sub>12</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = ε<sup>0</sup>x<sub>1</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = ε<sup>0</sup>x<sub>2</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>12</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>13</sub></td><td>ε<sup>0</sup><sub>13</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = ε<sup>0</sup>x<sub>1</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = ε<sup>0</sup>x<sub>3</sub>, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>13</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>23</sub></td><td>ε<sup>0</sup><sub>23</sub> = ε<sup>0</sup></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = ε<sup>0</sup>x<sub>2</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = ε<sup>0</sup>x<sub>3</sub>, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>23</sub>(ε<sup>0</sup>)<sup>2</sup></td></tr><tr><td>C<sub>44</sub></td><td>ε<sup>0</sup><sub>23</sub> = <sup>γ0</sup>&frasl;<sub>2</sub></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>3</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>2</sub>, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>23</sub>(<sup>γ0</sup>&frasl;<sub>2</sub>)<sup>2</sup></td></tr><tr><td>C<sub>55</sub></td><td>ε<sup>0</sup><sub>13</sub> = <sup>γ0</sup>&frasl;<sub>2</sub></td><td>u<sub>1</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>3</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>1</sub>, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>13</sub>(<sup>γ0</sup>&frasl;<sub>2</sub>)<sup>2</sup></td></tr><tr><td>C<sub>66</sub></td><td>ε<sup>0</sup><sub>12</sub> = <sup>γ0</sup>&frasl;<sub>2</sub></td><td>u<sub>1</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>2</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = <sup>γ0</sup>&frasl;<sub>2</sub>x<sub>1</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0</td><td>U<sub>e</sub> = <sup>V</sup>&frasl;<sub>2</sub>C<sub>12</sub>(<sup>γ0</sup>&frasl;<sub>2</sub>)<sup>2</sup></td></tr></tbody></table>

## 5.2. Stacking effect of GnR

The effective properties of the PVDF/GnR nanocomposites presented above were predicted under ideal conditions, in other words when all nano-shaped inclusions are well dispersed and aligned. However, the synthesis of nanocomposites is a challenging task and the ideal conditions studied above may be difficult to achieve. For example, graphene inclusions tend to stack up and to form graphite layers, leading to deterioration of the effective properties of the nanocomposites. The theoretical framework proposed in [30] was used to examine the effect of stacked GnR inclusions on the effective electro-elastic properties of nanocomposites. In the configuration presented in figure 7, all GnR inclusions are perfectly aligned with the third material direction of the PVDF matrix (i.e. the poling direction of PVDF). According to the literature, GnR inclusions can be considered as a quasi-1D nanomaterial. It was therefore
<table><caption>**Table 3.** Boundary conditions for piezoelectric constants [28]. The +and−respectively refer to the front face and back face of the RVE for a given specific direction.</caption><thead><tr><th>Piezoelectric constants</th><th>Deformation and applied electrical field</th><th>Displacement and applied electrical potentials</th><th>Electro-mechanical energy</th></tr></thead><tbody><tr><td>e<sub>31</sub></td><td>ε<sup>0</sup><sub>31</sub> = ε<sup>0</sup> ; E<sup>0</sup><sub>3</sub> = E<sup>0</sup></td><td>u<sub>1</sub>(+) = ε<sup>0</sup>x<sub>1</sub>, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0<br>ψ<sub>3</sub>(+) = 0, ψ<sub>3</sub>(−) = −E<sup>0</sup></td><td>U<sub>em</sub> = <sup>V</sup>&frasl;<sub>2</sub>e<sub>31</sub>ε<sup>0</sup>E<sup>0</sup></td></tr><tr><td>e<sub>32</sub></td><td>ε<sup>0</sup><sub>32</sub> = ε<sup>0</sup> ; E<sup>0</sup><sub>3</sub> = E<sup>0</sup></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = ε<sup>0</sup>x<sub>2</sub>, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = 0, u<sub>3</sub>(−) = 0<br>ψ<sub>3</sub>(+) = 0, ψ<sub>3</sub>(−) = −E<sup>0</sup></td><td>U<sub>em</sub> = <sup>V</sup>&frasl;<sub>2</sub>e<sub>32</sub>ε<sup>0</sup>E<sup>0</sup></td></tr><tr><td>e<sub>33</sub></td><td>ε<sup>0</sup><sub>33</sub> = ε<sup>0</sup> ; E<sup>0</sup><sub>3</sub> = E<sup>0</sup></td><td>u<sub>1</sub>(+) = 0, u<sub>1</sub>(−) = 0<br>u<sub>2</sub>(+) = 0, u<sub>2</sub>(−) = 0<br>u<sub>3</sub>(+) = ε<sup>0</sup>x<sub>3</sub>, u<sub>3</sub>(−) = 0<br>ψ<sub>3</sub>(+) = 0, ψ<sub>3</sub>(−) = −E<sup>0</sup></td><td>U<sub>em</sub> = <sup>V</sup>&frasl;<sub>2</sub>e<sub>33</sub>ε<sup>0</sup>E<sup>0</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Effective elastic properties predicted with MT and FE methods: Young's moduli (a) and shear moduli (b).</figcaption></figure>

assumed that the semi-major axis of the ellipsoids $a_3$ must tend towards infinity. The stacking direction is defined along the first direction of GnR inclusions and thus the aspect ratio $a_3/a_1$ of the ellipsoids is taken constant and equal to 10 000.0 (infinity). The second aspect ratio denoted by AR and defined by the semi-major axis ratio $a_2/a_1$ is used to depict the shape of the stacked GnR inclusions inside the nanocomposite. The AR values are ranged from 0.01 to 1.0; with a value of 0.01 theoretically corresponding to a single graphene layer for GnR inclusions. The numerical predictions presented in the section were obtained with MT micromechanics scheme.

As expected the inclusion shape has an effect on the Eshelby tensor, and hence the predicted values of effective elastic properties. In figure 7, the effective Young's modulus $Y_{11}$ values are slightly affected and a sharp decrease in the effective Young's modulus $Y_{22}$ values is predicted
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effective piezoelectric properties predicted with MT and FE methods.</figcaption></figure>

when the aspect ratio of GnR inclusions AR decreases to 0.1. For higher values of AR, the effective Young's moduli $Y_{11}$ and $Y_{22}$ are not sensitive to the values of the aspect ratio AR and to the GnR volume fraction, which is in agreement with the numerical study of Ji *et al* [30]. However, the effective Young's modulus $Y_{33}$ is not sensitive to the aspect ratio AR = $a_1/a_2$ and the so-called 'nano effect' seems to be preserved, figure 7(c). The predictions presented in figure 7 can be attributed to the particular shape of GnR inclusions (quasi 1D-geometry along the third direction). The values of the effective shear moduli $G_{13}$ and $G_{12}$ are not presented because they are nearly equal and are not altered by the stacking of GnR inclusions. In figure 8, the effective shear modulus $G_{23}$ strongly depends on the volume fraction of GnR and the stacking effect, i.e. the aspect ratio AR. However, when the aspect ratio of GnR becomes greater than 0.1, the effective shear modulus $G_{23}$ tends to a constant value independently of the GnR volume fraction values. In conclusion, the effective elastic constants are not altered by the addition of GnR for AR values greater than 0.1, with the exception of the effective Young's modulus $Y_{33}$.

The effective piezoelectric constants $e_{31}$ and $e_{33}$ predicted by the MT micromechanics method have fairly constant values regardless of the aspect ratio AR (i.e. the staking effect) and GnR volume fraction. The predicted values of the effective piezoelectric constant $e_{32}$ are presented in figure 9. The effective piezoelectric constant $e_{32}$ decreases as the aspect ratio AR of GnR inclusions tends to 0.1. Nevertheless, the stacking effect on the effective piezoelectric constants remains weak, which corroborates our assumption that GnR inclusions do not modify the $\beta$-phase fraction of PVDF matrix.

From the numerical results presented above, it can be concluded that GnR inclusions with an aspect ratio AR greater than 0.1 embedded in the PVDF matrix lead to a nanocomposite with the weakest mechanical behaviour. We have considered that a certain amount of GnR inclusions are stacked but have kept their initial orientation inside the PVDF matrix. The configuration is presented in figure 10 and the corresponding problem can be treated as a two-phase inclusion problem. The volume fraction of stacked GnR inclusions is denoted by $c_S$ and the total volume fraction of GnR inclusions by $c_R$. Following the method proposed in [30], the effective Young's moduli and piezoelectric constants of the nanocomposites with restacked
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Stacking effect of GnR inclusions on the normalized effective elastic properties: Young's modulus $Y_{11}$ (a), Young's modulus $Y_{22}$ (b) and Young's modulus $Y_{33}$ (c).</figcaption></figure>

GnR inclusions has then been computed with the micromechanics MT scheme. The predictions of the effective Young's modulus $Y_{22}$ and piezoelectric constant $e_{32}$ of the nanocomposites with restacked GnR inclusions are presented in figures 11 and 12 respectively.

In figure 11, the normalized values of the effective Young's modulus $Y_{22}$ decrease linearly as the volume fraction of stacked GnR inclusions (denoted $c_8$) increases. They also decrease as the GnR volume fraction decreases. The numerical predictions presented in figure 12 are in agreement with the predictions presented in [30]. However, no effect of stacked GnR is obtained for the two other effective Young's moduli, $Y_{11}$ and $Y_{33}$. As for the effective Young's modulus, the effects of stacked GnR inclusions are similar on the effective piezoelectric constants. In figure 12, the normalized values of the effective piezoelectric constant $e_{32}$ decrease linearly as the volume fraction of stacked GnR inclusions (denoted $c_5$) increases. The two other effective piezoelectric constants, $e_{31}$ and $e_{33}$, are almost independent of the volume fraction
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Stacking effect of GnR inclusions on the normalized effective shear modulus $G_{23}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Stacking effect of GnR inclusions on the normalized effective piezoelectric constant $e_{32}$.</figcaption></figure>

of stacked GnR inclusions and are not presented here. With the proposed micromechanical modeling, GnR stacking slightly causes to deteriorate the electro-mechanical behavior of the PVDF/GnR nanocomposite.

### 5.3. Alignment effect of GnR

The effects of piezoelectric inclusion orientation in infinite matrix have been well studied in previous numerical works on piezoelectric composites [25, 26]. For the different inclusions shapes studied in [26], the inclusions orientation has a significant effect on the effective electro-elastic behavior of the composites. In addition, the alignment effect of nano-inclusions on the electro-mechanical behavior of polymer nanocomposites was studied in [56] with a 2D-model
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Geometric configuration used in the micromechanics MT scheme with a two-phase inclusion problem to study the stacking effect of GnR.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Normalized effective Young's modulus $Y_{22}$ of GnR/PVDF nanocomposites where part of GNR inclusions are stacked.</figcaption></figure>

for the case of CNT nanocomposites. In the numerical predictions of [56], the electro-elastic behavior of the nanocomposite is improved when the CNT are all aligned along the strain loading direction. Since GnR can be considered as quasi 1D-geometry inclusions, it is expected that the electro-elastic behavior of GnR/PVDF nanocomposites is improved when all GnR are well aligned along the poling direction. Due to the anisotropy of the PVDF matrix and for sake of simplicity, only three geometric configurations were studied considering small rotations along the three material directions of the PVDF matrix. The various studied configurations are presented in figure 13. With the MT modeling previously presented, the alignment effect on the effective properties is only obtained in the configuration presented in figure 13(c).

The GnR alignment degree was expressed by the angle between the poling direction of the PVDF matrix (3-direction or third material direction) and the longitudinal axis of the GnR.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Normalized effective piezoelectric constant $e_{32}$ of GnR/PVDF nanocomposites where part of GnR inclusions are stacked.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** The different configurations studied: initial configuration of GnR within PVDF matrix (a), small rotation along $x_1$-direction (b), small rotation along $x_2$-direction (c) and small rotation along $x_3$-direction (d).</figcaption></figure>

The effective Young's moduli $Y_{11}$ and $Y_{22}$ are not influenced by small variations in the GnR alignment with respect to the three material directions of PVDF. This is the reason why numerical predictions for both effective Young's moduli are not presented. Because of the specific
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Alignment variations effect of GnR inclusions on the normalized effective Young's modulus $Y_{33}$.</figcaption></figure>

geometry of GnR (quasi-1D inclusion geometry), the effective Young's modulus $Y_{33}$ strongly depends on alignment variations by rotation along the 2-direction of PVDF matrix. The effective values of Young's modulus $Y_{33}$ are shown in figure 14. The effective Young's modulus $Y_{33}$ predicted for an alignment variation beyond 15° weakly depends on the value of the volume fraction of GnR in the PVDF matrix, figure 14. However, from the numerical predictions presented in figure 14, it can be concluded that the Young's modulus $Y_{33}$ can be improved significantly if the GnR inclusions length is aligned with an admissible variation of 5° along the material 3-direction of PVDF matrix (i.e. the poling direction of PVDF).

The effects of GnR alignment variations on piezoelectric constants $e_{31}$ and $e_{33}$ are quite negligible and are therefore not presented here. In figure 15, the normalized values of the effective piezoelectric constant $e_{32}$ decrease non-linearly as the alignment defect of GnR inclusions increases. However, the effective values of the piezoelectric constant $e_{32}$ remain close to zero. Therefore, it is not possible to conclude on the effects of the GnR alignment variations on the piezoelectric properties. This could be due to the absence of a model that describes the percolation behavior of GnR in the micromechanical modeling.

**6. Conclusion**

The effective elastic and piezoelectric properties of GnR/PVDF nanocomposites were studied using a micro-mechanics approach and taking into account the anisotropy of the PVDF matrix. The effective properties estimated with a micro-mechanics approach (MT scheme) were compared with the predictions provided by Finite Element simulations. For ideally aligned inclusions, a good agreement between the MT predictions and Finite Element predictions was obtained. From this preliminary study with perfectly aligned GnR inclusions, it was decided to study the effects of stacking and alignment of GnR inclusions on the effective properties of the nanocomposites. The stacking effect of GnR was modeled via the aspect ratio variations of the ellipsoidal inclusion in the MT micro-mechanics approach. In agreement with the literature, the stacking effect of GnR was found to have a great influence on the effective elastic
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Alignment variations effect of GnR inclusions on the normalized effective piezoelectric constant $e_{32}$.</figcaption></figure>

properties of the nanocomposite as it significantly improves the Young's moduli. However, the piezoelectric properties seem to be weakly impacted. The effects of varying GnR alignment were also studied using the MT approach. The improvement of elastic properties is strongly affected but the piezoelectric properties are not very dependent on GnR alignment variations. Our study leads to considering GnR as inclusions for designing PVDF nanocomposites with improved mechanical properties and uncompromised piezoelectric behavior.

**Data availability statement**

No new data were created or analysed in this study.

**Acknowledgment**

The authors gratefully acknowledge the funding from both parties: The French ANR ('Agence Nationale de Recherche' project EMBIA—'Electrical Multiscale Bioinspired Interfaces in Aeronautics', https://anr.fr/Project-ANR-21-CE05-0006) and the European Union's Horizon Europe research and innovation program (HORIZON-EIC-2023-PATHFINDEROPEN-01 under grant agreement No 101129952).

**Author contributions**

Salah Elbarnaty ID 0009-0008-8616-6381

Data curation (lead), Investigation (lead), Software (lead), Writing – original draft (lead)

Wiyao Azoti

Software (supporting), Writing – review & editing (supporting)

Joao Pedro M Correia ID 0000-0002-4881-8392

Conceptualization (equal), Supervision (lead), Writing – review & editing (lead)
Nadia Bahlouli

Funding acquisition (lead)

Saïd Ahzi

Resources (supporting), Supervision (supporting)

## Appendix A

Non-zero components of the Eshelby tensor S for the GnR-reinforced PVDF composite ($a_1 = 1$; $a_2 = 100$; $a_3 = 10000$) presented in the full 9 × 9 matrix form:

$$ S = \begin{bmatrix} 1.00126 & 0.739918 & 0.829358 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0.00226 & 0.006681 & 0.001934 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.001174 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.49866 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0.496495 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0.99681 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.002861 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} $$
Appendix B

```python
1 from odbAccess import *
2
3 # Open the ODB file
4 f = openOdb('C11_e31.odb')
5
6 # --- Downloading quantities needed for post-treatment ---
7
8 # Volume is the volume of each integration element of the RVE
9 Volume = f steps['Step-1'] .frames[-1].fieldOutputs['IVOL'] .values
10
11 # Epot is potential of each node of the RVE
12 EPot = f .steps['Step-1'] .frames[-1].fieldOutputs['EPOT'] .values
13
14 # ELSE is the elastic energy stored in each element of the RVE
15 ELSE = f .steps['Step-1'] .frames[-1].fieldOutputs['ELSE'] .values
16
17 # EENER is the dielectric energy stored in each element of the RVE
18 EENER = f .steps['Step-1'] .frames[-1].fieldOutputs['EENER'] .values
19
20 # --- Computing the total volume of the RVE ---
21 V = 0.0
22 for i in range(len(Volume)):
23     V += Volume[i].data
24
25 # --- Computing the sum of the Elastic and Dielectric energies ---
26 EE = 0.0  # Total elastic energy
27 EV = 0.0  # Total dielectric energy
28 for i in range(len(Volume)):
29     EE += ELSE[i].data
30     EV += EENER[i].data
31
32 # Em: total elastic energy stored in RVE
33 Em = EE * V
34
35 # Ed: total dielectric energy stored in RVE
36 Ed = EV * V
37
38 # Eem: total electroelastic energy stored in RVE
39 Eem = Em + Ed
40
41 # --- Parameters ---
42 eps_0 = 1e-01      # Strain used to create elastic energy
43 Voltage = 1e05     # Electric field between two opposed surfaces [V/m]
44
45 # --- Equations from Eq. (22),(23) ---
46 e31 = (2 * Eem) / (V * eps_0 * Voltage)
47 C11 = (2 * Em) / (V * (eps_0 ** 2))
48
49 # --- Display results ---
50 print('e31 =', e31)
51 print('C11 =', C11)
```

<caption>**Figure B1.** A typical python code to compute the elastic and piezoelectric coefficients $C_{11}$ and $e_{31}$ respectively on ABAQUS.</caption>
## References

[1] Maillard T, Claeyssen F, LeLetty R, Sosnicki O, Pages A and Vazquez Carazo A 2009 Piezomechatronic-based systems in aircraft, space, and defense applications *Proc. SPIE* **7331 73310K**

[2] Li M, Chen W M and Guan D 2004 Improvement of aircraft rolling power by use of piezoelectric actuators *Chin. J. Aeronaut.* **17** 87–92

[3] Zhu J, Yang J, Zhang W, Gu X and Zhou H 2023 Design and applications of morphing aircraft and their structures *Front. Mech. Eng.* **18** 34

[4] Marouf A, Charbonnier D, Vos J B, El Akoury R, Hoarau Y and Braza M 2024 Morphing of high-lift wing-flap system with cambering and trailing-edge flapping at high Reynolds number towards a full airplane application *J. Fluids Struct.* **127** 104111

[5] Abou-Khalil J et al 2024 Aerodynamic performance increase over an A320 morphing wing in transonic regime by numerical simulation at high Reynolds number *Int. J. Numer. Methods Heat Fluid Flow* **34** 10

[6] Krishnaswamy J A, Buroni F C, Garcia-Sanchez F, Melnik R, Rodriguez-Tembleque L and Saez A 2019 Improving the performance of lead-free piezoelectric composites by using polycrystalline inclusions and tuning the dielectric matrix environment *Smart Mater. Struct.* **28** 075032

[7] Krishnaswamy J A, Buroni F C, Garciá-Maciás E, Melnik R, Rodriguez-Tembleque L and Saez A 2020 Design of lead-free PVDF/CNT/BaTiO₃ piezocomposites for sensing and energy harvesting: the role of polycrystallinity, nanoadditives, and anisotropy *Smart Mater. Struct.* **29** 015021

[8] Kalimuldina G et al 2020 A review of piezoelectric PVDF film by electrospinning and its applications *Sensors* **20** 5214

[9] Sencadas V, Gregorio R and Lanceros-Méndez S 2009 α to β phase transformation and microstructural changes of PVDF films induced by uniaxial stretch *J. Macromol. Sci.* **48** 514–25

[10] Defebvin J 2015 Étude des relations structure-propriétés de matériaux hybrides piézoélectriques à base PVDF *PhD Thesis* University of Lille 1, France

[11] Talbourdet A 2018 Structures textiles piézoélectriques à base de PVDF pour la conversion d'énergie mécanique en énergie électrique *Textiles PhD Thesis* University of Lille, France

[12] Chiu F C, Chuang Y C, Liao S J and Chang Y H 2018 Comparison of PVDF/PVAc/GNP and PVDF/PVAc/CNT ternary nanocomposites: enhanced thermal/electrical properties and rigidity *Polym. Test* **65** 197–205

[13] Sodagar S, Saleh B, Fakhri P, Kashfi M, Feizi Mohazzab B and Momeni A 2020 Flexible piezoelectric PVDF/NDs nanocomposite films: improved electroactive properties at low concentration of nanofiller and numerical simulation using finite element method *J. Polym. Res.* **27** 203

[14] Xie L, Wang G, Jiang C, Yu F and Zhao X 2021 Properties and applications of flexible poly(Vinylidene fluoride)-based piezoelectric materials *Crystals* **11** 644

[15] Aqeel S M, Huang Z, Walton J, Baker C, Falkner D, Liu Z and Wang Z 2018 Polyvinylidene fluoride (PVDF)/polyacrylonitrile (PAN)/carbon nanotube nanocomposites for energy storage and conversion *Adv. Compos. Hybrid Mater.* **1** 185–92

[16] Chatterjee S, Nüesch F A and Chu B T T 2011 Comparing carbon nanotubes and graphene nanoplatelets as reinforcements in polyamide 12 composites *Nanotechnology* **22** 275714

[17] Wu C M, Chou M H and Zeng W Y 2018 Piezoelectric response of aligned electrospun polyvinylidene fluoride/carbon nanotube nanofibrous membranes *Nanomaterials* **8** 420

[18] Costa P, Silva J and Lanceros Mendez S 2016 Strong increase of the dielectric response of carbon nanotube/poly(vinylidene fluoride) composites induced by carbon nanotube type and pretreatment *Composites B* **93** 310–6

[19] Photopoulos P and Triantis D 2015 Dc conductivity measurements on PVDF composite samples of low grapheme content

[20] Ong M T and Reed E J 2012 Engineered piezoelectricity in graphene *ACS Nano* **6** 1387–94

[21] Chandratre S and Sharma P 2012 Coaxing graphene to be piezoelectric *Appl. Phys. Lett.* **100** 023114

[22] Young R J, Kinloch I A, Gong L and Novoselov K S 2012 The mechanics of graphene nanocomposites: a review *Compos. Sci. Technol.* **72** 1459–76

[23] Dunn M L and Taya M 1993 Micromechanics predictions of the effective electro elastic moduli of piezoelectric composites

[24] Dunn M L 1994 Electroelastic green's functions for transversely isotropic piezoelectric media and their application to the solution of inclusion and inhomogeneity problems *Int. J. Eng. Sci.* **32** 119–31
[25] Kuo W-S and Huang J 1997 On the effective electroelastic properties of piezoelectric composites containing spatially oriented inclusions *Int. J. Solid Struct.* **34** 2445–61

[26] Fakri N, Azrar L and El Bakkali L 2003 Electroelastic behavior modeling of piezoelectric composite materials containing spatially oriented reinforcements *Int. J. Solids Struct.* **40** 361–84

[27] Elouafi J, Azrar L and Aljinaidi A A 2015 Closed-form expressions for the effective moduli of heterogeneous piezoelectric materials *Int. J. Solids Struct.* **52** 19–32

[28] Odegard G M 2004 Constitutive modeling of piezoelectric polymer composites *Acta Mater.* **52** 5315–30

[29] Odegard G M, Gates T S, Wise K E, Park C and Siochi E J 2003 Constitutive modeling of nanotube-reinforced polymer composites *Compos. Sci. Technol.* **63** 1671–87

[30] Ji X Y, Cao Y P and Feng X Q 2010 Micromechanics prediction of the effective elastic moduli of graphene sheet-reinforced polymer nanocomposites *Modelling Simul. Mater. Sci. Eng.* **18** 045005

[31] Azrar L, Bakkali A and Aljinaidi A A 2014 Frequency and time viscoelectroelastic effective properties modeling of heterogeneous and multi-coated piezoelectric composite materials *Compos. Struct.* **113** 281–97

[32] Li J Y and Dunn M L 2001 Viscoelectroelastic behavior of heterogeneous piezoelectric solids *J. Appl. Phys.* **89** 2893–903

[33] Li K, Gao X L and Roy A K 2006 Micromechanical modeling of viscoelastic properties of carbon nanotube-reinforced polymer composites *Mech. Adv. Mater. Struct.* **13** 317–28

[34] Qin Q H 2004 Material properties of piezoelectric composites by BEM and homogenization method *Compos. Struct.* **66** 295–9

[35] Chwal M and Muc A 2021 FEM micromechanical modeling of nanocomposites with carbon nanotubes *Rev. Adv. Mater. Sci.* **60** 342–51

[36] Guo Z, Song L, Chai G B, Li Z, Li Y and Wang Z 2019 Multiscale finite element analyses on mechanical properties of graphene-reinforced composites *Mech. Adv. Mater. Struct.* **26** 1735–42

[37] Dutra T A, Ferreira R T L, Resende H B, Guimarães A and Guedes J M 2020 A complete implementation methodology for Asymptotic Homogenization using a finite element commercial software: preprocessing and postprocessing *Compos. Struct.* **245** 112305

[38] Chwal M and Muc A 2016 Transversely isotropic properties of carbon nanotube/polymer composites *Composites B* **88** 295–300

[39] Kumar D and Srivastava A 2016 Elastic properties of CNT-and graphene-reinforced nanocomposites using RVE *Steel Compos. Struct.* **21** 1085–103

[40] El Bahi A, Rouway M, Tarfaoui M, Moumen A E, Chakhchaoui N, Cherkaoui O and Omari L E H 2023 Mechanical homogenization of transversely isotropic CNT/GNP reinforced biocomposite for wind turbine blades: numerical and analytical study *J. Compos. Sci.* **7** 29

[41] Shingare K B and Naskar S 2021 Probing the prediction of effective properties for composite materials *Eur. J. Mech. A* **87** 104228

[42] Shang X, Liu Z, Zhang J, Lyu T and Zou Y 2023 Tailoring the mechanical properties of 3D microstructures: a deep learning and genetic algorithm inverse optimization framework *Mater. Today* **70** 71–81

[43] Del Castillo M and Pérez N 2021 Machine learning identification of piezoelectric properties *Materials* **14** 2405

[44] Wu J, Jiang J, Chen Q, Chatzigeorgiou G and Meraghni F 2023 Deep homogenization networks for elastic heterogeneous materials with two- and three-dimensional periodicity *Int. J. Solids Struct.* **284** 112521

[45] Liu R, Kumar A, Chen Z, Agrawal A, Sundararaghavan V and Choudhary A 2015 A predictive machine learning approach for microstructure optimization and materials design *Sci. Rep.* **5** 11551

[46] Singh K, Adhikari J and Roscow J 2024 Prediction of the electromechanical properties of a piezoelectric composite material through the artificial neural network *Mater. Today Commun.* **38** 108288

[47] Georgantzinos S K, Giannopoulos G I and Anifantis N K 2010 Numerical investigation of elastic mechanical properties of graphene structures *Mater. Des.* **31** 4646–54

[48] Vinogradov A and Holloway F 1999 Electro-mechanical properties of the piezoelectric polymer PVDF *Ferroelectrics* **226** 169–81

[49] Vinogradov A M and Holloway F 2000 Dynamic mechanical testing of the creep and relaxation properties of polyvinylidene fluoride *Polym. Test.* **19** 131–42
[50] Vinogradov A M, Schmidt V H, Tuthill G F and Bohannan G W 2004 Damping and electromechanical energy losses in the piezoelectric polymer PVDF *Mech. Mater.* **36** 1007–16

[51] Fan P, Wang L, Yang J, Chen F and Zhong M 2012 Graphene/poly(vinylidene fluoride) composites with high dielectric constant and low percolation threshold *Nanotechnology* **23** 365702

[52] Mano J F, Sencadas V, Costa A M and Lanceros-Méndez S 2004 Dynamic mechanical analysis and creep behaviour of β-PVDF films *Mater. Sci. Eng.* **370** 336–40

[53] Fang J, Vandenberghe W G and Fischetti M V 2016 Microscopic dielectric permittivities of graphene nanoribbons and graphene *Phys. Rev. B* **94**

[54] Eshelby J D 1957 The determination of the elastic field of an ellipsoidal inclusion, and related problems The determination of the elastic field of an ellipsoidal inclusion, and related problems *Proc. R. Soc. A* **241** 376–96

[55] Huang J H and Yu J S 1994 Electroelastic Eshelby tensors for an ellipsoidal piezoelectric inclusion *Compos. Eng.* **4** 1169–82

[56] Lee B M, Huang Z and Loh K J 2020 Effect of carbon nanotube alignment on nanocomposite sensing performance *Mater. Res. Express* **7** 046406