11851


FACTA UNIVERSITATIS  

Series: Mechanical Engineering  

https://doi.org/10.22190/10.22190/FUME050123059M 

© 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND 

Original scientific paper 

GEOMETRICALLY NONLINEAR ANALYSIS OF 

PIEZOELECTRIC ACTIVE LAMINATED SHELLS BY MEANS 

OF ISOGEOMETRIC FE FORMULATION 

Predrag Milić1, Dragan Marinković1,2, Žarko Ćojbašić1 

1University of Niš, Faculty of Mechanical Engineering in Niš, Serbia 
2Department of Structural Analysis, TU Berlin, Germany 

Abstract. The topic of piezoelectric active thin-walled structures has attracted a great 

deal of attention over the previous couple of decades. Lightweight structures with 

piezoelectric material based active elements, sensors and actuators, offer numerous 

advantages over their passive counterparts. This explains the motivation of authors to 

dedicate their work to this enticing research field. Accurate and reliable numerical 

tools for modeling and simulation of this type of structures is still a hot topic in the 

research community. This paper offers an isogeometric finite element formulation for 

shell type of structures made of composite laminates including piezoelectric layers 

characterized by the electro-mechanical coupling. The shell kinematics is based on the 

Mindlin-Reissner assumptions, thus including the transverse shear effects. A few 

examples selected from the available literature are considered to demonstrate the 

applicability of the developed numerical tool and assess its performance.   

Key words: Isogeometric analysis, Laminated structure, Reissner-Mindlin kinematics, 

Shell, Piezoelectricity, Geometrically nonlinear analysis 

1. INTRODUCTION 

Modern structures are made so as to be fit for their purpose. In order to meet this 

objective, they are designed to be lightweight. Modern fiber-reinforced composite 

materials enable design and manufacturing of lightweight structures, while improving the 

structural performance in different areas – the structures are characterized by high 

stiffness-to-weight ratio and become less expensive to operate with. Additionally, further 

improvement of structural performance is provided by incorporating active elements 

made of sophisticated multifunctional materials such as piezoelectric materials [1], shape 

                                                           
Received: January 05, 2023 / Accepted April 02, 2023 

Corresponding author: Predrag Milić  

University of Niš, Faculty of Mechanical Engineering in Niš, Aleksandra Medvedeva 14, 18000 Niš, Serbia,  

E-mail: predrag.milic@masfak.ni.ac.rs 



2 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

memory alloys [2], etc. The purpose of adding the active elements is to make the 

structures adaptive, and such structures are also often popularly referred to as “smart”. In 

this way, it is aimed at an essential change in their behavior from being only a passive 

player in the deformational behavior to being an active player that adapts to the current 

conditions. The structures based on this new paradigm include elements for monitoring 

structural parameters (sensors), elements that process signals (controllers), and elements 

that act onto the structure so as to influence their behavior [3]. The application of this 

type of structures is quite broad, thus covering various fields, including energy harvesting 

[4], robotics [5], space technologies [6], automotive [7], etc.  

A design based on thin walls is the most natural way of producing a lightweight 

structure, particularly if a structure is made of a fiber-reinforced composite. This 

approach, on the other hand, may increase lateral flexibility, thus significantly affecting 

deformational behavior of those structures, especially their dynamic behavior being 

affected [8]. Hence, the most common objective of adding the active elements is to 

positively influence the dynamic behavior of these structures with the aim of diminishing 

their vibrations and thus improve their safety and robustness, reduce noise emission, 

improve their life span, etc.  

The finite element method is currently one of the most advanced methods in the field 

of structural analysis. In the beginning, the developed elements were aimed at purely 

mechanical problems. In recent decades, multidisciplinary approaches gained in 

significance and the developments needed to cover coupled-field problems. The previous 

two decades have seen numerous developments of robust, accurate and efficient 

formulations and finite elements for modeling piezoelectric structures characterized by 

electromechanical coupled field effects [9, 10, 11, 12].  

Various 2D theories have been developed and applied to modeling and simulation of 

thin-walled structures [13]. While the classical laminate theory is appealing due to less 

nodal degrees of freedom required, the higher order of continuity of the shape functions is 

the major obstacle in its implementation into the FEM. Despite of that, the issue was 

addressed in a discretized form by some authors [14]. Most of the authors turned their 

attention to the first-order shear deformation theory, which is due to both, the better 

suitability of the theory for laminated structures (prone to the transverse shear effects) and 

the required C0-continuity of the shape functions [15, 16]. Taking an improved accuracy 

and a more consistent description of the transverse shear effects as the main objectives, 

some authors based their developments on higher-order shear deformation theories [17, 

18]. The numerical effort associated with these theories, which is not far away from the 

full 3D approach, is the obvious disadvantage. Besides these efforts, the need to cover 

geometrically nonlinear effects was recognized recently and a number of developments 

were reported in this direction [19, 20]. Recent papers were dedicated to the co-rotational 

FE formulation for both passive and active shell structures [21, 22].  

Another relatively recent and quite interesting approach is denoted as isogeometric 

analysis. It is originally proposed by Hughes et al. [23] and aims at seamless integration 

of the CAD geometry into the FE models. The most distinctive property of this approach 

is that the actual, i.e. CAD geometry remains unchanged upon the FE discretization. 

Development of isogeometric FE formulations relying on both the classical laminate 

theory based on the Kirchhoff-Love kinematics [24, 25], and the first-order shear 

deformation theory based on the Mindlin-Reissner kinematics [26], were reported. The 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 3 

same is valid for higher-order theories [27]. Isogeometric FE formulations for composite 

laminates were also reported [28, 29]. 

This paper aims at a geometrically nonlinear isogeometric formulation for shell type 

of structures made of composite laminates with embedded piezoelectric layers that may 

act both as actuators and sensors. It represents an extension of the recently reported 

isogeometric development for the same type of structures, but done for linear analysis 

[30].  

2. NURBS GEOMETRY, BASIC FUNCTIONS AND MESH PROPERTIES 

CAD description of geometry and, therewith, the isogeometric formulation of finite 

elements is based on Non-Uniform Rational B-Splines – NURBS. NURBS is a special 

type of B-spline and, hence, its properties partly depend on the properties of B-spline. To 

determine the basis functions of the NURBS, it is necessary to define the basis functions 

of the B-spline first. A clear advantage of NURBS compared to a B-spline resides in its 

ability of accurate representation of complex geometric shapes.  

The basis functions of B-spline of the order zero are determined as: 

  
1

,0

1

0 otherwise

i i

i
N


 

 


  
  (1) 

A B-spline of the order p is defined based on the Cox-De Boor's recursive formula. 

For this purpose, a knot vector =[0, 1,..., n+p+1] is needed. Higher order basic functions 
are defined as follows: 

      
1

, , 1 1, 1

1 1

i pi
i p i p i p

i p i i p i

N N N
  

  
   

 

  

   


 

 
 (2) 

A set of parametric coordinates i given in a non-decreasing order defines the knot 
vector. The continuity of the basis functions is of great importance and it is defined by the 

function order, p, and the knot multiplicity, k. A knot multiplicity equals to k implies Cp-k 

continuity. In addition to this, Ni,0() is a stepped function giving non-zero values only in 

a half-open interval [i, i+1). Furthermore, Ni,p() is given by a linear combination of 
two functions of degree (p-1). A basic function of the order p gives non-zero values only 

in a half-open interval [i, i+p+1) and the sum of all basic functions of the order p at 

any point  is equal to 1. Finally, the basis functions are linearly independent and non-
negative. 

A NURBS curve of the order p is defined by using the B-spline basis functions and 

the control polygon points with the corresponding weight coefficients or with control 

polygon points and their corresponding NURBS basis functions: 

 
,

0
,

0
,

0

( )

( ) ( )

( )

n

i p i i n
i

i p in
i

i p i
i

N w P

C R P a b

N w



  









   





 (3) 



4 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

Here, Pi are the control polygon points, wi the corresponding weights, {Ni,p()} the p-

order B-spline basic functions defined on the non-uniform knot vector ={a,...,a, p+1,..., 

m-p-1, b,...,b}, while Ri,p represent the basic rational functions of the order p that form 
NURBS. 

Typically, the knot vector is normalized. This implies that the first index coordinate 

of the vector is 0 (a = 0), while the last one is 1 (b = 1). The number of repetitions of the 

elements a and b in the knot vector defines the order of spline. Using this parameter, one 

cab realize the discontinuity at the spline ends. Hence, the following equation defines a 

NURBS surface of the orders p and q in the - and -directions, respectively: 

 

, ,
0 0

0 0
, , ,

0 0

( ) ( )

( , ) ( , ) 0 , 1

( ) ( )

n m

i p j q ij ij n m
i j

ij ijn m
i j

k p l q k l
k l

N N w P

S R P

N N w

 

     

 

 

 

 

   





 (4) 

where Rij are the rational B-spline basis functions of NURBS surface, i.e. the NURBS 

basis functions: 

 
, ,

, ,
0 0

( ) ( )
( , ) 0 , 1

( ) ( )

i p j q ij

ij n m

k p l q kl
k l

N N w
R

N N w

 
   

 
 

  


 (5) 

Fig. 1 shows an example of a NURBS surface (patch) given by Eq. (4). Points of the 

control polygon in the physical space, index vectors in the index space and the basic 

functions order define the shown NURBS surface. The element boundaries are set in the 

parameter space. The surface belonging to an element is defined on half-open intervals 

[i, i+1) and [j, j+1) and they are recognizable after the mapping to the physical 
space (Fig. 1). 

 

Fig. 1 NURBS surface – mapping between the physical and parametric space 

The surface of an element is only a part of the whole patch. It is defined by the 

necessary mesh parameters, which include control polygon points, weight coefficients and 

the NURBS basis functions belonging to the half-open intervals [i, i+1) and [j, 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 5 

j+1). The NURBS basic functions belonging to the points of the patch control polygon 
outside the observed half-open interval yield a zero value. Consequently, they have no 

influence on this part of the geometry. In order to keep track of this fact, all the control 

points that define the element geometry and their basic functions are numbered from 1 to 

nen, with nen=(p+1)(q+1).  

The initial mesh of finite elements has a total number of finite elements determined 

by the number of elements in the index vectors and the degrees of the basic functions. It 

also implies that there is no repetition of elements in the index vectors except at the 

beginning and the end of the vector: 

   en n p m q    (6) 

Here, m and n are the number of points of the control polygon along the - and - 

direction, respectively, while p and q are the basic function degrees in the - and - 
directions, respectively. There are several possibilities to change the initial mesh. Knots 

may be inserted or removed from the knot vectors (h-refinement), the degree of basic 

functions may be changed (p-refinement) and the degree of basic functions may also be 

changed followed by insertion of new knots into the knot vectors (k-refinement). The 

degree of continuity of the basic functions across the elements may be affected by the last 

mentioned technique.  

3. ISOGEOMETRIC SHELL FORMULATION 

Four coordinate systems are used in the formulation of the shell isogeometric finite 

element. Those would be the global Cartesian coordinate system (x, y, z), the natural 

coordinate system (r, s, t) with coordinate values -1<r, s, t<+1, the local Cartesian 

coordinate system (x’, y’, z’), and the curvilinear coordinate system (, , ) (Fig. 2). 

The tangent plane at any point of the surface is given by vectors 
1t
'V and 

2t
'

'
V  

determined as derivatives of the position vector r with respect to and

 

'
1

'
2

t
1

'

t
1

( , )

( , )

en

en

n
k

k
k

n
k

k
k

Rr
V P

Rr
V P

 

 

 

 






 
 


 
 





 (7) 

  

Fig. 2 Coordinate systems used in shell element formulation 



6 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

Here, Rk(,) are the basic functions and Pk the control polygon points on the element 

surface. The resulting vectors, 
1t
'V and

2t
'

'
V , are not perpendicular to each other. In order to 

generate an orthogonal coordinate system, the vector product of the obtained two vectors 

is computed to give a vector normal to the tangent plane.  

 ' ' '
1 2

'

t tz
V V V   (8) 

The unity vectors of the coordinate system are then easily determined:  

 
2 1t t
' ' '

z
V V V   (9) 

 
' ' '
1 2

'' '' ''

' ' '
1 2

1 2 3

t t

1 2 3

t t
1 2 3

= , = , = ,

'' '' ''

z'' '' ''

x y z
'' '' ''

z

l l l
V V V

e m e m e m
V V V

n n n

     
     

       
     
     

 (10) 

While the definition of the normal to the mid-surface is a trivial task in classical 

FEM, the isogeometric formulation is somewhat challenging with respect to this aspect. 

The challenge originates from the fact that the control polygon points are not placed on 

the actual shell surface and they don’t have a unique projection onto the shell surface. The 

issue has been addressed in the recent work by Milic et al. [30], where several methods 

from the literature are reported and the approach based on the Greville points [31] is 

adopted. The same approach is also used in this work.  

A point on the structure that is not on the mid-surface is determined by the position 

vector of its projection onto the mid- surface and with the distance in the direction 

perpendicular to the mid-surface. The thickness of the shell is determined interpolating 

the known thicknesses hk at the polygon control points by means of their corresponding 

NURBS functions. The position of any point is given as: 

 

3

3
=1

3

( , )
2

     
      

      
     
      


en

k k xn
k

k k k y
k

k k z

x x V
h

y R y t V

z z V

   (11) 

Here, xk, yk, zk represent the control point global coordinates, Vk3x, Vk3y, Vk3z are the 

components of the vector normal to the mid-surface at control point Pk.  

The same shape functions are applied to interpolate the displacement field and define 

the element geometry. Hence: 

 
=1

en

R

k kn
R

k k k
k R

k k

u u u

v R v v

w w w

    
     

      
     
      

  (12) 

Here,  uk, vk and wk are the control point displacements in the x-, y- and z directions 

respectively, while uRk, v
R

k and w
R

k denote the relative displacements caused by rotations 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 7 

at the control point k. The latter are to be expressed via the nodal rotations xk, yk and zk, 
which represent the global rotations. 

The Reissner-Mindlin kinematics is applied to determine the relative displacement of 

any point within the shell with respect to its counterpart on the mid-surface, Those are 

defined with respect to the local coordinate system at control point k, as follows: 

 
l

2 l 2 l

2
2 2

   
        

             
        

  

R

k xk

x'kkR k k
k k k yk k k

y'kR k
k zk

u
Vh h

v t -V ,V t -V ,V
V

w









 (13) 

Here x’k and y’k are the nodal rotations with respect to the x'- and y'-axis Fig. 4. 
Material properties of composite laminates (Fig. 3) are directionally dependent. This 

fact calls for the definition of the strain field with respect to the local coordinate system. 

This is commonly done so that a distinction is made between the in-plane components 

(membrane-flexural strains) and the out-of-plane components (transverse shear), as given 

below: 

  

 

 

'

ε

ε

'

x'x'

' y'y'
mf

x'y'

'

s
y'z'

x'z

u'

x'

v'

y'

u' v'

y' x'

v' w'

z' y'

u' w'

z' x'











 
 


 

 
 

  
                
          

       
      

    
      

 
  

 
  

 (14) 

  

Fig. 3 Laminated structure 

The derivatives of the local displacements with respect to the local coordinates are 

obtained from the derivatives of the global displacements in the global coordinate system 

by means of the transformation matrix:  



8 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

 

l

2 l 2 3

3

, , ,

, ,

, , ,

, , ,

k,x ,y ,x ,x ,x ,x

, ' ,

,y' k ,y ,y ,y k k k,y ,y

, , ,
,z ,z ,zk,z ,z ,z

u v w V u v w

u v w V u v w V V V

u v wVu v w

     
     

         
     

      

 (17) 

where, for instance, ,
,

,x
u = ∂u/∂x and ,xu  = ∂u/∂x. 

The derivatives are transformed from the natural to the global coordinate system 

using the Jacobian inverse matrix: 

  
, , ,

, , ,
=1

, , ,

3 3 3
2 2 2

en

k k k
k k k

n
k k k

k k k
k

t t t

k k k
k k x k k y k k z

N N N
x y z

x y z
N N N

J x y z x y z

x y z
h h h

N V N V N V

  

  

  

  

   
 

  
  
    

    
  

  
   

 
 

  (18) 

  
, , ,

1

, , ,

,x ,x ,x

,y ,y ,y

,z ,z ,z ,t ,t ,t

u v w u v w

u v w J u v w

u v w u v w

  

  



   
   

   
   
   

 (19) 

Upon all the mentioned transformations, the strain field given in Eq. (16) has the 

following form: 

      
 

 0 1

mf Tm Rf
T

u

s Ts R s R s R

B B | t B u

' = B u = - - u - - + - - - -

B B | B + t B u



                 
       
    
                     

 (20) 

The element strain–displacement matrix, [Bu] consists of the membrane-flexural, [Bmf], 

and transverse shear, [Bs], strain–displacement matrices, which further involve 

translations and rotations related parts (denoted by T and R in the subscript).  

4. CONSTITUTIVE RELATIONS OF PIEZOLAYERS 

The form of piezoelectric constitutive equations depends on the choice of 

independent variables. Within the framework of FEM, the strain and the electric filed are 

the usual choice, thus yielding the matrix form of piezoelectric constitutive equations: 

 
       

      

T
-

-

E

e

= C e E

D = e d E

 



 
 

 
 

 (21) 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 9 

Here {σ} and {ε} are the mechanical stress and strain in the Voigt notation, respectively, 

[CE ] is the symmetric material constitutive (Hook’s) matrix, {D} and {E} are the 

dielectric displacement and the electric field vector, respectively, [e] is the piezoelectric 

coupling matrix, and [dε] is the vector of dielectric constants.  

For quite thin piezopatches, it has been shown that the electric field may be assumed 

to be constant over the thickness of the piezolayers [32]. Hence, for the kth piezolayer: 

 
k

k

k

E
h


  (22) 

where k is the difference of electric potentials between the electrodes of the k
th 

piezolayer and hk is the thickness of the piezolayer. The resulting  electric field – electric 

potential matrix has a quite simple, diagonal form: 

 
1

k

B
h



 
 
 
 
 
  

 (23) 

5. GEOMETRICALLY  NONLINEAR FINITE ELEMENT ANALYSIS 

Linear structural analysis resides on the assumption that the displacements are 

relatively small compared to the dimensions of the modeled structure, while the material 

behavior is also linear. Also, the boundary conditions must remain unchanged over the 

course of deformation. On the other hand, geometrically nonlinear analysis takes into 

account the change in structural configuration during the deformation. Deformation of the 

structure leads to a change in structural parameters. Considering the fact that that the 

structural deformation progresses continuously, the structural parameters also change are 

continuously. In each structural configuration, it is necessary to establish the equilibrium 

between the generalized internal and external forces.  

The solution of a nonlinear analysis is sought in a step-by-step approach. If the 

equilibrium of the system at time t is known, a solution at time t+t is sought with a 

suitably chosen time step t. Even in a static analysis, a label t is established as an 

auxiliary variable that indicates the progress of the load level. A set of linearized FE 

equations for the piezoelectric continuum are given on the element level as follows: 

 
       

       

T e e

e e

t t t tt tt
euu e u ext int

t t t tt tt
eu e ext int

K u K f f

K u K q q



 









         

          

 (24) 

Here, {fexte}, {finte}, {qexte}and {qinte} are the external and internal mechanical and 

electrical loads, [KuuT], [Ku] and [K] are the tangential mechanical stiffness, 

piezoelectric coupling and dielectric stiffness matrices, respectively, while the vectors 

{ue} and { e } include the mechanical and electrical degrees of freedom of the element, 



10 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

respectively. All the terms are defined for time t, which is denoted by the left superscript, 

and  is used to denote the increment of the corresponding quantity between two 

configurations. The external excitations are the applied forces/moments and electrical 

charges. Their internal counterparts are obtained by integrating the mechanical stresses 

and dielectric displacements over the current configuration of the structure, respectively: 

      
T

t
e

t t t

int LV
f B dV   (25) 

    
T

t
e

t t t

int A
q B D dA     (26) 

The integration in Eq. (25) runs over the whole structure, V, while in Eq. (26) only over 

the surface of the piezopatches, A. 

Updated Lagrangian formulation is used in thgis work and, hence, all the vectors and 

matrices use the last determined configuration as a reference configuration. The tangential 

stiffness matrix is computed using the following equation: 

      
T T

t t

t t t t t tE

uuT u uV V
K B C B dV G G dV                (27) 

where t[G] is the matrix of partial derivatives and t[] is the Cauchy stress tensor: 

    
     
     
     

T 0 0

0 0

0 0

g

g G g
x y z

g

 
     

     
     

 

 (28) 

  

   

   

   

0 0

0 0

0 0

xx xy xz

xy yy yz

xz yz zz


  

     

  


  
    

              
   
      

 (29) 

The piezoelectric stiffness matrix, [Ku], and the dielectric stiffness matrix, [K], are 

computed in the same manner as in the linear analysis, with the only difference that they 

are integration over the current volume of the structure: 

 
T

t

t t t

V
K B d B dV


  

               (30) 

  
TT

t

t t tt

u u uV
K B e B dV K

  
                (31) 

It should be emphasized that the piezo-electrically induced loads are configuration 

dependent. They are obtained by multiplying the piezoelectric coupling matrix, Eq. (31), 

with the electric potential.  

The obtained nonlinear problem is resolved by means of the Newton–Raphson 

iterative procedure [33]. 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 11 

6. NUMERICAL EXAMPLES 

In this section, a few selected numerical examples of laminar piezo-electric shells 

under different load sets are considered. Bothe the actuator and sensor role of 

piezoelectric layers are covered.  

6.1 Nonlinear analysis of a plate 

The in-plane dimension of considered plate are 4040 mm, while the thickness varies 

in the considered cases. The structure was discretized using the following finite elements: 

Shell91 (Ansys), ACShell9 [34] and the developed formulation (NURBS). In its initial 

geometry, the structure is flat, but becomes curved after the first loading step. First, the 

results of linear and geometrically nonlinear analysis of the plate made of only one 

piezoelectric ceramic isotropic layer PZT G1195 (properties in Table 1) with different 

types of finite elements are compared. The structure is exposed to an external force of 24 

kN in all considered cases. In the first two cases, the force is distributed over a 40 mm 

long line (600 kN/m) and in the third case, it is a concentrated force (Fig. 5).  

Table 1 Material properties [34] 

Material 

properties 

PZT G1195 

piezolayer 

T300/976 

graphite/epoxy 

Elastic properties 

Y11 (GPa) 63 150 

Y22 (GPa) 63 9 

12 0.3 0.3 

G12 (GPa) 24.2 7.1 

Piezoelectric properties 

e31 (10-5 C/mm2) 2.286 0 

e32 (10-5 C/mm2) 2.286 0 

 

 

Fig. 5 Initial flat structures, boundary conditions and loads 

Figure 5 shows the three models of the considered plate. Due to the kinematic and 

dynamic boundary conditions, the maximal displacement in the first case would be larger 

than in the latter two for the same material properties and the thickness. The first model is 

a clamped plate with a thickness of 5 mm with a load at the free end. The thickness of the 

plate is chosen so that an obvious non-linear effect occurs. The results of linear and non-

linear analysis of this model are given in Fig. 6. Fig. 6a shows the vertical displacement of 

the free edge mid-point (point A), while Fig. 6b gives its displacement in the x-direction. 

The latter result is equal to zero in the linear analysis. 



12 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

 

Fig. 6 Clamped structure-displacement of free edge mid-point in: a) vertical direction, 

and b) in longitudinal (x-)direction 

The latter two cases deal with the structure simply supported over the opposite two 

edges and with the structure simply supported over all four edges. The position of the 

load was chosen adequately so that, in the first case, it is a line distributed force acting 

over the middle span width and, in the second case, it is a concentrated force acting at the 

structure mid-point. Due to the boundary conditions, the structural behavior is 

significantly stiffer compared to the first case. In order to have geometrically nonlinear 

effects, a smaller thickness of 3 mm is chosen. In the second example, point B, i.e. the 

free edge mid-point is chosen as the point with the largest transverse displacement. The 

maximal transverse displacement of point B in linear and nonlinear analysis is given in 

Fig. 7a. The linear result for the displacement in x- and y-direction of any point of the 

plate is zero. However, the nonlinear analysis gives a non-zero displacement in the y-

direction, as shown in Fig. 7b. 

 

Fig. 7 2-edge-simply supported structure displacement of free edge mid-point in: 

a) vertical direction and b) in longitudinal (y-)direction 

In the case of a structure with four simply supported edges, the mid-point of the plate 

(point C) was chosen as the reference, Fig. 5c. Due to the susceptibility to shear locking 

effect of the Ansys Shell91 element, the initial mesh (4x4) was corrected to 2020 

elements to provide accuracy comparable with the ACShell9 and NURBS finite elements, 

which use a 44 mesh. Both, linear and nonlinear prediction yield only the transverse 

deflection of the observed point Fig 8.  



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 13 

 

Fig. 8 4-edges-simply supported structure: structure mid-point vertical displacement 

In all three considered cases the initial geometry is a plate and the applied loads are 

such that the linear prediction involves only flexural and shearing strains (no membrane 

strains). However, with progressing nonlinearity, the membrane strains are also induced 

as the configuration changes are taken into account. It is obvious and expected that the 

effect of nonlinearity is the least pronounced in the case of a clamped structure. The 

model with the NURBS elements in the previously shown examples has a very good 

agreement with the results obtained with the other considered elements. 

The first two plate models are further used to illustrate piezoelectric actuation of the 

structure in the realm of nonlinear deformations. The observed structures are considered 

to be made of three layers. The outer layers are made of PZT material (Tab. 1) and the 

middle layer is made of graphite-epoxy material with a fiber orientation of 90 with 

respect to the global x-axis. Given that the piezo layer has a relatively small actuation, the 

aforementioned models should be modified in order to produce a geometrically nonlinear 

effect. For this reason, a thickness of 0.15 mm is assigned to the piezo layers have and the 

thickness of 0.2 mm to the composite layer. The piezolayers are oppositely polarized, so 

that their actuation forms a moment equally distributed along the edges of the plate. A 

voltage of 300V is supplied to both piezolayers. The fiber orientation is chosen so that the 

structure has a lower bending stiffness about the y-axis than the bending stiffness about 

the x-axis.  

The induced piezo-electric coupling forces are configuration dependent. The new 

increment of the actuating bending moments is calculated based on the piezo-electric 

coupling stiffness matrix by performing its integration over the actual structural 

configuration. The direction and intensity of the induced loads depends on the current 

configuration (follower forces). The follower nature of the induced forces/moments 

demands small increments [33] of the electric voltage. 

An approximate approach is also used by treating the cases as purely mechanical with 

the loads that are not of the follower type. To achieve this, a full voltage of 300 V is 

applied to the initial structure and the bending moments are computed. The calculated 

moment is then applied in increments.The calculations done with the ACShell9 and 

NURBS finite element yield small differences (in the order of 10-1 %). 



14 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

 

Fig. 9 Initial geometry and boundary conditions  

The diagrams in Fig. 10 are obtained by means of the ACShell9 and NURBS diagrams 

using incremental voltages. The purely mechanical case with he predefined moments is 

coputed using the Shell91 element. Those results are also included in the diagrams. In the 

clamped structure model (Fig. 9a), two points, A and B, on the structure are chosen to 

monitor their motion with respect to the gradually increasing excitation. Figs. 10a and 10b 

depict the transverse deflection and the displacement in the x-direction of the point A. 

Fig. 11 shows equivalent results to those given in Fig. 10, but for point B of the structure. 

 

Fig. 10 Clamped piezoelectric plate, point A displacements: a) y-direction; b) x-direction  

 

Fig. 11 Clamped piezoelectric plate, point B displacements: a) y-direction; b) x-direction 

The previously discussed example reveals a small nonlinear effect on the 

displacements in the transverse deflection. Fig. 12 gives the results for the transverse 



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 15 

deflection in the case shown in Fig. 9b. Obviously, the difference between the linear and 

nonlinear results is significantly larger in this case, which is the consequence of more 

pronounced nonlinear effect through the induced membrane strains and stresses. 

 

Fig. 12 2-edges-simply supported piezoelectric structure: displacement of point C in the 

vertical direction 

Observing the results of the examples shown, a negligible difference can be observed 

between the application of elements in which the voltage was applied incrementally 

(ACShell9, NURBS) and elements, in which the precomputed moment was applied 

incrementally (Ansys Shell91). This should not lead to a general conclusion in the same 

direction.  

6.2 Nonlinear analysis of piezo-laminated semicircular arch 

In the next example, a composite cylinder arc was considered. The structure is formed 

from three layers. The middle passive layer is made of metal with a thickness of 5.842 

mm (Y=68.95 GPa, =0.3, =7750kg/m3). The outer layers are piezoelectric layers with a 

thickness of 0.254 mm (Y=63 GPa, =0.3, =7600 kg/m3, e31=e32=16.11 C/m
2, 

d33=1.6510
-8 F/m). 

 

Fig. 13 Piezo-laminated semicircular arch  



16 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

Only one force F, the magnitude of which is 100 N, acts onto the structure in the 

middle of the tip of the free edge. The response of the structure in the nonlinear static 

analysis under the action of the force F was observed through the displacement of the tip 

of the free edge in the radial direction as well as in the circular direction. The acting force 

also induces electric voltages in the piezoelectric layers. The voltage of the inner 

piezoelectric layer is observed for the comparison purposes with the results from the 

literature. Actually, Tzou and Ye [35] proposed the example originally. It was further 

modified by Zhang [36], and Marinkovic et al. [22]. Zhang discretized the structure with 

10 elements and using one element in the width direction. The model by Marinkovic at al   

[22] model uses 160 triangle elements. The NURBS model was formed with quadratic 

basic functions, with 4 elements in the width direction and 20 elements in the hoop 

direction. The displacements are shown in Fig. 14 and the induced electric potentials of 

the inner piezolayer are given in Fig. 15. 

 

Fig. 14 Displacement of the free arc edge tip in the radial and hoop directions  

 

Fig. 15 Sensor voltage of the inner piezolayer  



 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 17 

Minor differences between the results for the induced electric potential in the inner 

piezolayer can be noticed in Fig. 15. This can be attributed to different meshes applied 

and also the fact that Zhang [36] used the rigorous geometrically nonlinear formulation, 

while Marinkovic et al. [22] used a co-rotational formulation. 

7. CONCLUSIONS 

Shell type of structures, i.e. thin-walled structures represent a natural way of achieving 

lightweight design. In addition to this, novel structural materials, such as fiber-reinforced 

laminates with high stiffness-to-weight ratio, are used to produce even lighter structures 

with an option of tailoring material properties by a careful selection of the composite 

material composition, number and sequence of layers as well as the fiber orientations in 

the layers. This has been the state-of-the-art in many high performance structures for a 

few decades. A further improvement comes in the form of adaptive structures – structures 

provided with active elements, i.e. sensors and actuators, and which can proactively react 

to different conditions they are exposed to. In the specific case addressed in this paper, 

the active elements are in the form of thin piezoelectric layers embedded in the laminate 

structure.  

This paper proposed a novel numerical tool for modeling the mentioned type of 

structures. The finite element method is addressed as the most powerful method in the 

field of structural analysis. The development is based on a relatively novel isogeometric 

approach. The main advantage that it offers is a seamless integration of the CAD 

geometry into the FE model, thus keeping the original, correct geometry of the structure. 

The proposed formulation covers the mechanical field in the whole structure as well as 

the electrical field and the piezoelectric coupling in the piezoelectric layers. Thin-walled 

structures are prone to deformations characterized by large local rotations whereby the 

strains remain small. For that reason, the developed tool covers both linear and 

geometrically nonlinear analysis. For the mechanical field, an equivalent single layer 

approach based on the First-order Shear Deformation Theory is applied.  

The considered cases are selected from the available literature in order to have 

reference solutions. They involve flat and curved structures. But it should be emphasized 

that although the considered structure is flat in the first set of examples, it is essentially 

dealt with a shell structure in the geometrically nonlinear analysis as the deformed 

structure is curved. The computations were done for both purely mechanical and electric 

excitations. Also, the actuator and sensor function of active layers were covered. It was 

shown that the proposed formulation successfully covers all the mentioned features and 

represents a good match for the elements based on the classical FE formulation.  

In the future work, the performance of the developed isogeometric formulation with 

the mesh refinement should be addressed. Also, the formulation should be extended to 

cover material nonlinearities in the piezoelectric layers that occur at higher voltages.  

Acknowledgement: This research was financially supported by the Ministry of Science, 

Technological Development and Innovation of the Republic of Serbia (Contract No. 451-03-

47/2023-01/ 200109) and by the Science Fund of the Republic of Serbia (Serbian Science and 

Diaspora Collaboration Program, Grant No. 6497585).  



18 P. MILIĆ, D. MARINKOVIĆ, Ž. ĆOJBAŠIĆ 

REFERENCES 

1. Wang, K., Alaluf, D., Rodrigues, G., Preumont, A., 2021, Precision Shape Control of Ultra-thin Shells with 

Strain Actuators, Journal of Applied and Computational Mechanics, 7(Special Issue), pp. 1130-1137. 

2. Todorov, T., Mitrev, R., Penev, I., 2020, Force analysis and kinematic optimization of a fluid valve driven by 

shape memory alloys, Reports in Mechanical Engineering, 1(1), pp. 61-76. 

3. Preumont, A., 2018, Vibration Control of Active Structures, Springer Cham, 518 p. 

4. Mitrev, R., Todorov, T., Fursov, A., Fomichev, V., Il'in, A., 2021, A Case Study of Combined Application of 

Smart Materials in a Thermal Energy Harvester with Vibrating Action, Journal of Applied and Computational 

Mechanics, 7(1), pp. 372-381. 

5. Nguyen, X., Nguyen, H., 2022, Investigation of influences of fabrication tolerances on operational 

characteristics of piezo-actuated stick-slip micro-drives, Facta Universitatis-Series Mechanical Engineering, 

20(1), pp. 109-126. 

6. Shao, S., Song, S., Xu, M., Jiang, W., 2018, Mechanically reconfigurable reflector for future smart space 

antenna application, Smart Materials and Structures, 27(9), 095014. 

7. Kulkarni, H., Zohaib, K., Khusru, A., Shravan-Aiyappa, K., 2018, Application of piezoelectric technology in 

automotive systems, Materials Today: Proceedings, 5(10, Part 1), pp. 21299-21304. 

8. Nestorović, T., Marinković, D., Chandrashekar, G., Marinković, Z., Trajkov, M., 2012, Implementation of a user 

defined piezoelectric shell element for analysis of active structures, Finite Elements in Analysis and Design, 52, 

pp. 11-22. 

9. Rama, G., Marinković, D., Zehn, M., 2018, Efficient three-node finite shell element for linear and geometrically 

nonlinear analyses of piezoelectric laminated structures, Journal of Intelligent Material Systems and Structures, 

29(3), pp, 345-357. 

10. Carrera, E., Valvano, S., Kulikov, G.M., 2018, Electro-mechanical analysis of composite and sandwich 

multilayered structures by shell elements with node-dependent kinematics, International Journal of Smart and 

Nano Materials, 9(1), pp. 1-33. 

11. Rama, G., Marinkovic, D.Z., Zehn, M.W., 2017, Linear shell elements for active piezoelectric laminates, Smart 

Structures and Systems, 20(6), pp. 729-737. 

12. Kulikov, G.M., Plotnikova, S.V., Carrera, E., 2018, Hybrid-Mixed Solid-Shell Element for Stress Analysis of 

Laminated Piezoelectric Shells through Higher-Order Theories, in: Altenbach, H., Carrera, E., Kulikov, G. 

(Eds), Analysis and Modelling of Advanced Structures and Smart Systems, Advanced Structured Materials, vol 

81. Springer, Singapore. 

13. Reddy, J.N., 2003, Mechanics of Laminated Composite Plates and Shells: Theory and Analysis, Second 

Edition. CRC Press, 858 p. 

14. Gabbert, U., Koppe, H., Seeger, F., Berger, H., 2002, Modelling of smart composite shell structures, Journal of 

Theoretical and Applied Mechanics, 3(40), pp. 575-593. 

15. Rama, G., Marinkovic, D., Zehn, M., 2018, A three-node shell element based on the discrete shear gap and 

assumed natural deviatoric strain approaches, Journal of the Brazilian Society of Mechanical Sciences and 

Engineering, 40(7), 356. 

16. Petrolo, M., Carrera, E., 2021, Selection of element-wise shell kinematics using neural networks, Computers & 

Structures, 244, 106425. 

17. Carrera, E., Zozulya, V.V., 2022, Carrera unified formulation (CUF) for the micropolar plates and shells. I. 

Higher order theory, Mechanics of Advanced Materials and Structures, 29(6), pp. 773-795.  

18. Li, G., Carrera, E., Hou, Y., Kulikov, G.M., 2021, Multi-layered plate finite element models with node-

dependent kinematics for smart structures with piezoelectric components, Chinese Journal of Aeronautics, 

34(8), pp. 164-175. 

19. Rama, G., Marinkovic, D., 2018, High performance 3-node shell element for linear and geometrically 

nonlinear analysis of composite laminates, Composites Part B: Engineering, 151, pp. 118-126. 

20. Katariya, P.V., Hirwani, C.K., Panda, S.K., 2019, Geometrically nonlinear deflection and stress analysis of 

skew sandwich shell panel using higher-order theory, Engineering with Computers, 35, pp. 467–485. 

21. Liguori, F.S., Madeo, A., 2021, A corotational mixed flat shell finite element for the efficient geometrically 

nonlinear analysis of laminated composite structures, International Journal for Numerical Methods in 

Engineering, 122(17), pp. 4575-4608. 

22. Marinković, D., Rama, G., Zehn, M., 2019, Abaqus implementation of a corotational piezoelectric 3-node shell 

element with drilling degree of freedom, Facta Universitatis-Series Mechanical Engineering, 17(2), pp. 269- 

283. 

https://www.tandfonline.com/journals/umcm20
https://www.sciencedirect.com/journal/chinese-journal-of-aeronautics
https://www.sciencedirect.com/journal/composites-part-b-engineering
https://link.springer.com/journal/366


 Geometrically Nonlinear Analysis of Piezoelectric Active Laminated Shells by Means of Isogeometric... 19 

23. Hughes, T.J.R., Cottrell, J.A., Bazilevs, Y., 2005, Isogeometric analysis: CAD, finite elements, NURBS, exact 

geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering, 194(39-41), pp. 

4135- 4195. 

24. Nguyen-Thanha, N., Valizadehb, N., Nguyenc, M.N., Nguyen-Xuand, H., Zhuange, X., Areiasf, P., Zih, G., 

Bazilevsg, Y., De Lorenzisa, L., Rabczukb, T., 2015, An extended isogeometric thin shell analysis based on 

Kirchhoff-Love theory, Computer Methods in Applied Mechanics and Engineering, 284, pp. 265-291. 

25. Milić, P., Marinković, D., 2015, Isogeometric FE analysis of complex thin-walled structures, Transactions of 

FAMENA, 39(1), pp. 15-26. 

26. Benson, D.J., Bazilevs, Y., Hsu, M.C., Hughes, T.J.R., 2010, Isogeometric shell analysis: The Reissner-Mindlin 

shell, Computer Methods in Applied Mechanics and Engineering, 199(5-8), pp. 276-289. 

27. Echter, R., Oesterle, B., Bischoff, M., 2013, A hierarchic family of isogeometric shell finite elements, Computer 

Methods in Applied Mechanics and Engineering, 254, pp. 170-180. 

28. Yujie, G., Martin, R., 2015, A layerwise isogeometric approach for NURBS-derived laminate composite shells, 

Composite Structures, 124, pp. 300-309.  

29. Hosseini, S., Remmers, J.J., Verhoosel, C.V., de Borst, R., 2015, Propagation of delamination in composite 

materials with isogeometric continuum shell elements, Int. J. Numer. Methods Eng., 102(3), pp. 159-179. 

30. Milić, P., Marinković, D., Klinge, S., Ćojbašić, Ž., 2023, Reissner-Mindlin based isogeometric finite element 

formulation for piezoelectric active laminated shells, Tehnički Vjesnik, 30(2), pp. 416-425.  

31. Adam, C., Bouabdallah, S., Zarroug, M., Maitournam, H., 2015, Improved numerical integration for locking 

treatment in isogeometric structural elements, Part II: Plates and shells, Computer Methods in Applied 

Mechanics and Engineering, 284, pp. 106-137. 

32. Marinković, D., Köppe, H., Gabbert, U., 2007, Accurate modeling of the electric field within 

piezoelectric layers for active composite structures. Journal of Intelligent Material Systems and 

Structures, 18(5), pp. 503-513. 

33. Bathe, K.J., 1996, Finite Element Procedures, Prentice Hall, New York. 

34. Marinković, D., Koppe, H., Gabbert, U., 2008, Degenerated shell element for geometrically nonlinear analysis 

of thin-walled piezoelectric active structures, Smart materials and structures, 17(1), 015030. 

35. Tzou, H.S., Ye, R., 1996, Analysis of piezoelastic structures with laminated piezoelectric triangle shell 

elements, AIAA J, 34(1), pp. 110–115. 

36. Zhang, S., 2014, Nonlinear FE Simulation and Active Vibration Control of Piezoelectric Laminated Thin-

Walled Smart Structures, PhD Thesis, Fakultat fur Maschinenwesen der Rheinisch-Westfalischen Technischen 

Hochschule Aachen, Germay, 210 p.