416 Technical Gazette 30, 2(2023), 416-425
ISSN 1330-3651 (Print), ISSN 1848-6339 (Online) https://doi.org/10.17559/TV-20230128000280
Original scientific paper
Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active
Laminated Shells
Predrag MILIĆ, Dragan MARINKOVIĆ*, Sandra KLINGE, Žarko ĆOJBAŠIĆ
Abstract: The paper deals with the isogeometric analysis (IGA) of active composite laminates with piezoelectric layers. IGA is a special formulation of the finite element
method (FEM) that aims at seamless integration of geometric and finite element modelling. NURBS basis functions are employed to develop isogeometric shell formulation
based on the Reissner-Mindlin kinematics. Piezolayers characterized by electro-mechanical coupled field effects enable active behavior of the considered structures. The
electric field acts across the thickness of the piezolayers and is coupled to the in-plane strains. In addition to a number of advantages that NURBS modelling provides,
defining the surface normal vector at the points of the control polygon, which are generally not located on the surface, creates certain difficulties. A method of determining
the surface normal vectors at the points of the control polygon based on the Greville's points is discussed. In order to demonstrate the applicability of the developed
formulation, a benchmark case is computed and the results are compared with those obtained by means of classical FEM formulation, which are available in the literature.
Keywords: grevill points; isogeometric analysis; laminated Reissner-Mindlin shell; piezolayer
1 INTRODUCTION
The progress of science and technology led to the
emergence of modern structures, which are not only
lightweight, but also kind of 'smart'. The lightweight
property is attributed to advanced materials, primarily
fiber-reinforced composites, while the structure 'smartness'
is due to the so-called active materials like piezoelectric
materials [1, 2], shape memory alloys [3], etc. These new
materials enabled a dramatic change in the behavior of the
structures from conventional passive ones to actively
controlled structures, which include capabilities such as
self-sensing, signal processing, and control by means of
actuation [4]. Smart structures are used in various fields of
engineering: aeronautics and space engineering [5], energy
harvesting [6, 7], medical applications [8], robotics [9], etc.
Particularly lightweight structures benefit from this new
design approach in order to reduce the mass of thin-walled
structures, the thickness is often reduced to the limits of
unwanted deformation, noise and vibrations, which calls
for the control of their dynamic behavior [10].
The focus of this paper is on modelling thin-walled
active structures with embedded piezoelectric layers,
which offer both sensing and actuating options by means
of electro-mechanical coupling. To investigate the
kinematics of laminated shell structures, various 2D
theories can be used [11]. Some developments were based
on the classical laminate theory [12]. However, the first-
order shear deformation theory was much more employed
due to its better suitability for laminated structures, which
are prone to transverse shear effects [13, 14]. Higher-order
shear deformation theories were applied as well [15, 16],
however the numerical effort with these theories is nearly
the same as with the full 3D approach, which is
disadvantageous. Novel developments were focused on
geometrically nonlinear formulations [17, 18].
Isogeometric analysis, proposed by Hughes et al. [19],
aims to close the gap between the CAD (or actual)
geometry and the geometry used upon finite element (FE)
discretization. There are several different isogeometric
shell formulations: Kirchhoff-Love [20-22], Mindlin-
Reissner [23], higher-order kinematics [24]. Models of
isogeometric laminated composite shells have also been
developed [25, 26].
The accuracy of the results of the isogeometric
analysis of the shell largely depends on the precisely
determined normal vectors at the points of the control
polygon. Considering that they do not have a unique
projection on the middle plane of the shell, several
researchers have investigated this problem [27, 28].
Upon providing the motivation to work on this topic
and a short overview of different directions of development
in the field of active piezoelectric laminates, the rest of the
manuscript is organized as follows. The next section
provides basics related to B-splines and NURBS basic
functions. This is followed by the procedure of finite
element mesh refinement and the determination of the
director vectors normal to the reference surface at control
points, both within the framework of isogeometric FE
formulation. Upon defining the displacement field, the
mechanical and piezoelectric constitutive equations are
used to define the finite element matrices for the coupled
field problem. The derived FE formulation is applied to a
benchmark case in order to demonstrate the development.
Finally, the concluding remarks of the study are given.
2 B-SPLINE AND NURBS BASIC FUNCTIONS
In this work, the geometry of the FE model is based on
NURBS surfaces (NURBS-Non-Uniform Rational
B-Splines). The model properties are influenced by this
fact and, therefore, the basis spline (B-spline), NURBS
functions, as well as curves and surfaces based on them are
addressed below.
NURBS has a clear advantage compared to a B-spline
related to the possibility of accurate descriptions of quite
complex geometric shapes. NURBS are actually built by
using B-splines as basic functions.
Cox-De Boor's recursive formula [29] is used to define
a B-spline of the order p based on a knot vector
= [
0,
1, ...,
n+p+1]. For p = 0, the basic functions read:
1
0
1
0otherwise
ii
i,
N
(1)
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
Tehnički vjesnik 30, 2(2023), 416-425 417
Higher order basic functions are defined in this
manner:
1
,,1 1,1
11
ip
i
ip ip i p
ip i ip i
NN N
(2)
The knot vector is given as a set of parametric
coordinates
i in a non-decreasing order. The function
order, p, and knot multiplicity determine the continuity of
basis functions. If the knot multiplicity is k, the continuity
is Cp−k. Further characteristics of the basic functions are:
Ni,0(
) is a stepped function equal to zero for all
except
for a half open interval
[
i,
i+1); Ni,p(
) is defined as a
linear combination of two functions of degree (p − 1); a
basic function of order p has a value different from zero
only in the semi-interval
[
i,
i+p+1); sum of all basic
functions of the order p at a point
is equal to 1 (partition
of unity); the functions are non-negative and linearly
independent.
A p-order NURBS curve is represented as a rational
function based B-spline basic functions as given below
[29, 30]:
,
0
,
0
,
0
,
n
ip i i n
i
ip i
n
i
ip i
i
NwP
CRPab
Nw
(3)
where Pi are the control points forming a control polygon,
wi are the weights, {Ni,p(
)} are the B-spline basic
functions of the order p defined on the non-uniform knot
vector = {a, ..., a,
p+1, ...,
m-p-1, b, ..., b}, and Ri,p are the
p-order basic rational functions of the NURBS.
The knot vector is usually normalized, so that a = 0 and
b = 1. The order of spline defines how many times the
elements a and b are repeated in the knot vector. In this
way, the discontinuity is realized at the ends of the spline.
A NURBS surface of the order p in the
direction and
order q in the
direction is given as:
,,
00
,,,
00
, 0 , 1
nm
i p j q ij ij
ij
nm
kp lq kl
kl
NNwP
S
NNw
(4)
Introducing the basic rational function:
,,
,,
00
, 0 , 1
ip jq ij
ij nm
kp lq kl
kl
NNw
R
NNw
(5)
the equation for the surface reads:
00
, , , 0 , 1
nm
ij ij
ij
SRP
(6)
3 FINITE ELEMENT MESH AND DIRECTOR VECTORS AT
CONTROL POINTS
An example of a NURBS surface (patch) defined as
per Eq. (5) is shown in Fig. 1. The NURBS surface is
defined by points of control polygon in the physical space,
index vectors in the index space and by degrees of the basic
functions. The boundaries of element surfaces are defined
in the parameter space (Fig. 2). Upon mapping into the
physical space (Fig. 1), the elements of the NURBS surface
can be observed. Each element-surface is defined on half
open intervals
[
i,
i+1); and
[
j,
j+1). This is the
initial mesh of elements.
Figure 1 NURBS surface based on quadratic basic functions
The geometry of the element reference surface (as part
of the patch) is determined by the mesh of the control
polygon points with the corresponding weight coefficients
and corresponding NURBS basis functions
,
ij
R
on
the half open intervals
[
i,
i+1) and
[
j,
j+1). Other
NURBS basic functions of the patch outside the half open
interval have a zero value, and therefore the points of the
control polygon corresponding to them have no influence
on this part of geometry. For the sake of simplicity, all
control points that affect the geometry of element as well
as their corresponding basic functions will be numbered
from 1 to nen, where nen = (p + 1)(q + 1). Fig. 1 shows a
patch with the element and its corresponding control
polygon points marked.
Figure 2 Index space, parametric space and natural coordinate system
Isogeometric analysis closes the gap between the CAD
geometry obtained from a CAD modeller and the finite
element geometry. The initial number of elements along
each direction is determined for the initial mesh by the
following equation and is clearly visible in the parameter
space:
e
nnpmq (7)
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
418 Technical Gazette 30, 2(2023), 416-425
where m and n are the number of basic functions (i.e. points
of control polygon) along the
and
directions, and p and
k are the degrees of the basic functions in the
and
directions.
The initial mesh can be changed with the techniques of
knot insertion into the knot vectors (h-refinement),
increasing the degree of basic functions (p-refinement) and
increasing the degree of basic functions with the
subsequent insertion of a new knot into the knot vectors
(k-refinement). With the last mentioned technique, the
degree of the basic functions and the density of the mesh
are increased with the aim of increasing the degree of
continuity at the elements boundaries.
For the formulation of the shell isogeometric finite
element type, we need four coordinate systems:
- the global Cartesian coordinate system (x, y, z),
- the natural coordinate system (r, s, t) - an orthonormal
system with −1 < r, s, t < +1,
- the local Cartesian coordinate system (x', y', z'),
- the curvilinear coordinate system (
,
,
).
Figure 3 Defining the coordinate system at a point on the reference surface
The tangent plane at a point on the surface is defined
by vectors
1
t'
V
and
2
t'
V
. They are determined as derivatives
of the position vector r by the variables
and
.
11
21
,
,
nen k
'
'k
tk
nen k
'
'k
tk
R
r
VP
R
r
VP
(8)
where Rk(
,
) and Pk are the basic functions and the points
of the control polygon on the element reference surface on
which the observed point is located. The obtained vectors,
1
'
Vt
and 2
'
Vt
, are not generally perpendicular to each other.
It is necessary to form an orthogonal coordinate system.
The vector product of the previous two vectors gives a
vector normal to the tangent plane.
12
'
'''
ztt
VVV
(9)
An orthonormalized coordinate system was formed by
Eq. (11):
21
'''
tzt
VVV
(10)
123
12
'' 1 '' 2 '' 3
123
12
=, =, =
'' '' ''
'''
ttz
'' '' ''
xyz
'' '' ''
'''
ttz
lll
VVV
ememe m
VVV
nnn
(11)
In order to form a finite element model with a shell
element type using the isogeometric approach, it is
necessary to determine the vectors normal to the reference
surface at the points of the control polygon. In the
isogeometric finite element method, most of the control
polygon points do not lie on the reference surface, in
contrast to the classical finite element method where the
nodes are located on the reference surface. In addition, they
do not have a unique projection on the reference surface.
To determine the vector normal to the reference surface for
a certain point of the control polygon, we can find several
methods in the literature:
- basis systems obtained by closest point projection,
- calculation of exact basis systems,
- determination of normal vectors in Greville points.
The first method of determining the base system is
based on determining the shortest distance between the
point of the control polygon and the reference surface. In
general, it is not possible to obtain a solution in closed form
for this problem. The iterative Newton-Raphson method is
most commonly used to solve it. Once the projection point
is determined, the orthonormalized coordinate system can
be determined according to the previous procedure.
The second method for the exact basis systems
calculation uses the fact that, for each element, the normal
vectors to the reference surface can be determined at the
integration points in two ways. First, the vector normal to
the reference surface at the integration point can be
determined through the derivative of the position vector
with the respect to the variable
or
(as shown earlier)
and, second, by interpolating over the normal vector at the
points of the control polygon (which are currently
unknown). It follows that a system of (p + 1)(q + 1)
equations has to be formed for each component of all
element control points. Also in this case the minimum
number of integration points is (p + 1)(k + 1). By increasing
the degree of the basic functions, the calculation of the
normal vector of the point of the control polygon to the
reference surface becomes more difficult. For example, for
the degree of basic functions 4, the number of points of the
control polygon that is active per element is 25. The
minimum number of integration points (to solve the system
of equations for the element) per direction is 5, and the total
number of equations would be 75 for all three components
of the normal vector. In addition, in this way it is necessary
to determine one more component of the basic system. If
the model is more complex, i.e. with more elements, the
problem becomes more demanding to solve.
The third method, which is used here, is based on the
determination of the control polygon normal vectors to the
reference surface using the Greville points. Greville points,
or abscissas, have one-to-one correspondence with the
control polygon points. Greville points are defined in the
parametric space as follows [31]:
++
GP GP
mmp nnq
mn
,
pq
(12)
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
Tehnički vjesnik 30, 2(2023), 416-425 419
The number of Greville points corresponds to the
number of NURBS basis functions. Hence, it is also equal
to the number of control points. Tangent vectors to the
reference surface can be formed at Greville's points:
1
2
GP GP
=1
GP GP
=1
,
,
en
'
en
'
nkm n
GP
k
t
k
nkm n
'GP
k
t
k
R*
r
VP
R
r
VP
(13)
In the previous equations, k denotes the number of the
element control point and basis function. Rk represents the
basic functions of the element that have a value different
from zero, Pk represents the points of the control polygon
of the corresponding basic functions Rk. The vectors
formed at the Greville points lie in the tangent plane to the
reference surface at the point with parametric coordinates
GP GP
,
mn
. Given that we know the parametric coordinates
and normal vectors in the Greville points, we can form a
system of equations based on interpolation. By solving the
system of equations, we determine the components of
vector normal to the reference surface corresponding to the
point of control polygon. The point of the control polygon
does not have a direct projection on the reference plane.
The reason for this is the property of the basic functions of
NRUBS, not all of which have the maximum value of 1
(only those on the patch boundaries have the maximum
value of 1). The components of the vector normal to the
reference surface in the Greville points can be calculated
based on the known other two vectors:
12
12
GP GP
GP
3GP GP
''
''
''
tt
''
tt
VV
V
VV
(14)
On the other hand, the vector components in the
Greville points can be represented by the following
equation:
'
GP GP GP
=1
, , 1, 2, 3
en
n
n
km n kl
t
k
VR Vl
(15)
In the previous equation, the vectors Vkl at the control
points are unknown. In practice, integration is often
performed with a reduced number of integration points. In
that case, it is not possible to determine the normal vectors,
so the determination of normals is done using Greville's
points.
4 DIPLACEMENT FIELD
The position of the point on the middle surface can be
presented as a function of parametric coordinates:
=1
(, )
en k
n
kk
k
k
x
x
y
Ry
zz
(16)
The thickness of the shell is assumed to be in the
direction normal to the reference surface of the element and
can be obtained for any point of the reference surface by
means of the interpolation functions and the thickness hk
corresponding to each control point.
3
3
=1
3
(, ) 2
en kkx
n
k
kk rkky
k
kkz
xxV
h
yR y ttV
zzV
(17)
where xk, yk, zk are the global coordinates of the control
point of element, hk is the control point correspondent
thickness; trk is the offset of the reference surface from the
mid-surface in the natural c.s. of control point k, given as
trk = 2hr/hk; Vk3x, Vk3y, Vk3z are component of the vector
perpendicular to the reference surface at control point Pk
(Eq. (15)).
The isogeometric formulation of the element implies
the same shape functions for the interpolation of the
displacements field as for the CAD and element geometry.
Thus, it will be:
=1
en
R
kk
n
R
kk k
kR
kk
uuu
vRvv
www
(18)
where uk, vk and wk represent the nodal displacements along
x, y and z axis (global displacements) respectively and uRk,
vRk and wRk are the relative displacements of the point on
the thickness direction line through the control point k with
respect to the control point, resulting from the rotation od
the line, and of course along x, y and z axis. The
displacements uRk, vRk and wRk are to be expressed in terms
of the nodal rotations
xk,
yk and
zk (global rotations).
Figure 4 Equivalent layer approach for multilayer material and coordinate systems
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
420 Technical Gazette 30, 2(2023), 416-425
In order to do that, the assumption that the thickness
direction line remains straight after deformation will be used
(Reissner-Mindlin kinematics). Threfore, in the local
coordinate system fixed to the control point k, the following
relations can be written:
,
2
R
kxk
kl
Rk
krkk2klyk
Rk2
kzk
uq
V
h
vt-t-VV q
V
wq
(19)
θ
θ
θ
xk
x'k kl
yk
y'k k2
zk
qV
qV
(20)
where x'k and y'k are the nodal local rotation (about the x'
and y' axis) Fig. 5.
Figure 5 Rotations of the thickness direction line in the local coordinate system
Directionally dependent material properties demand to
define the strain field in the local coordinate system. This
is usually done so that a distinction is made between the in
plane (membrane-flexural) components and the out-of-
plane (transverse shear) components:
'
ε
ε
'
x'x'
'y'y'
mf
x'y'
'
sy'z'
x'z
u'
x'
v'
y'
u' v'
y' x'
v' w'
z' y'
u' w'
z' x'
e
e
e
(21)
Derivatives of displacements in the natural coordinate
system can be represented by transformation using
derivatives of global displacements in the global
coordinate system and the transformation matrix:
l
2l23
3
'' '
,x' ,x' ,x'
'' '
,y' ,y' ,y'
'' '
,z' ,z' ,z'
k,x ,x ,x
k,y,y,ykkk
,z ,z ,z
k
uvw
uvw
uvw
Vuvw
VuvwVVV
uvw
V
(22)
where, for example, u', x' = ∂u'/∂x' and u, x = ∂u/∂x.
The transformation of derivatives from the natural to
the global coordinate system is done by means of Jacobian
inverse matrix:
,,,
,,,
,,,
=1
333
,, ,
1
,, ,
222
en
ttt
kkk
kkk
n
kkk
kkk
k
kkk
kkx kky kkz
,x ,x ,x
,y ,y ,y
,z ,z ,z ,t ,
xyz
xyz
xyz
NNN
xyz
NNN
xyz
hhh
NV NV NV
uvw uvw
uvw uvw
uvw uv
J
J
t,t
w
(23)
Using the local displacement derivatives given in Eq.
(22), the following form is obtained:
01
mf
u
s
Tm Rf T
Ts R s R s R
'= = --
|t
-- + -- --
|+t
B
eBu u
B
BB
u
BBBu
(24)
Obviously, the strain–displacement matrix, [Bu]
contains the membrane-flexural, [Bmf], and transverse
shear, [Bs], strain displacement matrices.
5 PIEZOELECTRIC CONSTITUTIVE RELATIONS
Selecting the strain and the electric field as
independent variables, the linear piezoelectric constitutive
equation in the matrix form is:
ε
T
E
=
=
CεeE
DeεdE
σ
(25)
where [CE] is the symmetric material constitutive (Hook's)
matrix, {
} is the stress vector, {
} is the strain vector, [e]
is the piezoelectric coupling matrix, {E} the electric field
vector, {D} the dielectric displacements vector, {d} the
vector of dielectric constants (electric displacement to
electric field ratio). Stress and strain are second order
tensors, but also that for the sake of convenience, they are
in the engineering (Voight) notation commonly
represented in the form of vectors {
} and {
}.
Material properties of commercially available
piezoelectric ceramics are transversely isotropic.
The above equations can be given together in a developed
form, as follows:
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
Tehnički vjesnik 30, 2(2023), 416-425 421
'
'
'
11 11 12 13 31
12 11 13 31
22
13 13 33 31
33
55 15
23
55 15
13
11 12
12
115
2
3
σ
σ
σ
σ
σ
1
σ2
''
''
''
''
''
''
CCC 0 0 0 |0 0 e
CCC 0 0 0 |0 0 e
CCC 0 0 0 |0 0 e
000C 0 0 |0 e 0
0000C 0 |e 0 0
00000 CC |0 0 0
----- - |- - -
D0000e 0 |d
D
D
11
22
33
23
13
12
1
11
15 11 2
31 31 33 33 3
ε
ε
ε
γ
γ
γ
''
''
''
''
''
''
'
'
'
E
00
E
000e 0 0 |0 d 0
eee 00 0 |0 0 d E
(26)
The electrodes are placed on the upper and lower
surface of the piezolayers. It is assumed that the electric
charge is uniformly distributed over the electrodes and,
consequently, that the electric field acts in the thickness
direction only, so that E1' = E2' = 0. The components of the
dielectric displacements D1' and D2' are not necessarily
equal to zero. However their product with the
corresponding components of the electric field (E1' and E2',
respectively) in the electric Gibbs energy is zero, which
allows the seventh and the eight row and column to be left
out of consideration. Hence, the constitutive equation of
the piezoelectric layer reads:
11 11
11 12 31
22 22
12 11 31
12 12
66
55
23 23
55
13 13
31 31 33
33
000
000
00 00 0
000 0 0
0000 0
000
'' ''
'' ''
'' ''
'' ''
'' ''
' '
QQ |e'
QQ |e'
Q|
Q|
Q|
|
e' e' | d'
DE
(27)
where the following reduced coefficients are introduced:
22
13 13
11 11 12 12
33 33
2
15
55 55 66 11 12
11
2
13 33
31 31 33 33 33
33 33
Q
1
Q2
'
CC
QC , C ,
CC
e
QC , CC
d
Ce
ee e,d'd ,
CC
(28)
It is easy to notice that all the shear strain and stress
components are decoupled from the rest in Eq. (27).
However, this is generally not the case with the in-plane
shear components for the considered passive material
(fiber-reinforced composite laminates). In order to keep the
same approach to obtaining the cross-sectional force and
moment resultants as for the passive material, the equations
for the transverse shear components will be separated from
the in-plane components, thus yielding:
11 11
11 12 31
22 22
12 11 31
66
12 1 2
31 31 33
33
σε
σε
σγ
'' ''
'' ''
'' ' '
''
QQ 0|e'
QQ 0|e'
00Q|0
---|-
e' e' 0 | d'
DE
(29)
23 23
13 13
σε
σγ
'' ''
'' ''
55
kQ
(30)
where k is the shear correction factor.
The electric field distribution in the thickness direction
of piezoelectric layers should satisfy both, the Maxwell's
equations for dielectrics and the element kinematics. The
Reissner-Mindlin kinematics would imply linear electric
field across the thickness [32]. However, the same
investigation [32] shows that the classical assumption of
constant electric field over the thickness offers satisfying
accuracy for typical, quite thin piezopatches. Hence, the
electric field of the kth piezolayer is given by:
k
k
k
Eh
(31)
with
k denoting the difference of electric potentials
between the electrodes of the kth piezolayer and hk is the
thickness of the piezolayer. This yields a diagonal form of
the electric field electric potential matrix:
1
k
Bh
(32)
6 FINITE ELEMENT DISCRETIZED EQUATIONS
Discretization of a continuum results in a finite
element representation. Proceeding in the way
characteristic for the finite element method, the generalized
displacement field that includes the mechanical
displacement field {u}T = {u, v, w}T in the global
coordinate system (x, y, z) and the electric potential
of
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
422 Technical Gazette 30, 2(2023), 416-425
each piezolayer are related to the corresponding nodal
values {u}i and i by means of interpolation functions.
In equilibrium, the virtual work of all acting forces in
moving through a virtual displacement is zero:
δ+δ0
IE
WW (33)
Assuming a system of conservative forces and that the
electric charge of the system is conserved as well
(electrically isolated system), the previous equation yields:
1
S1
T
TTT
Tε
TT
1
T
2
11
δδδ
d
δ
+δdδd
δδdδ0
S2
E
V
vS
V
nm
Pjj
jj
jj
V
VS
fq S Q ,
εCεεeE Eeε
EdE
uF uF
uF
(34)
where Fs1 are the external forces acting on surface S1, q is
the external charge defined on surface S2, and Qj represent
the electric concentrated charge at m points.
The element mechanical stiffness matrix resulting
from the variational principle reads [33]:
uu mf s t
=++
KKKK
(35)
The element stiffness matrix includes membrane-
flexural (in plane) [Kmf], transverse shear [Ks] and the so-
called torsional (drilling) stiffness [Kt]. The former two are
easily computed based on the derived strain field:
T
1e
eT
1
T
s
TT
01
e
T
01e
d
t
d
Tm
mf m Tm R f
V
Rf
Ts
s
Rs Rs
V
Ts R s R s
|t V
t
|tV
B
KCBB
B
B
C
KBB
BB B
(36)
where [Cm] is the part of Hook's matrix relating the in-plane
stresses and strains, and [Cs] relates the transverse shear
stress and strains.
A node has five degrees of freedom in the local
coordinate system (two rotations), but as a consequence of
transformation there are in general six degrees of freedom
in the global coordinate system. This may give rise to
numerical problems for a special position of the element
leading to zero stiffness for the rotation about one of the
global axes. An erratic behavior of the element.
A numerically quite simple solution is proposed by
Zienkiewicz and Taylor [34] and it implies introduction of
additional torsional stiffness. Governing torsional strain
energy, Et, is defined and it serves as a penalty function
forcing the local rotation z' to be approximately equal to
0.5(∂v'/∂x' − ∂u'/∂y') at integration points:
2
n
()
11''
θd
22''
n
tz'
Ar,s,o
vu
EYh A
xy
(37)
Here, Y is the Young's elasticity module and
n is a
fictitious parameter that can be chosen.
The element dielectric stiffness matrix is given as:
Td
V
V
KBdB (38)
It is diagonal as the difference of the electric potentials
of one piezoelectric layer affects only the electric field of
the very same piezolayer. Analytical integration in the
thickness direction yields:
11
33
1
11
1
0
det d d
2
0
'
k
''
kk
drs
hzz
KJ
(39)
where J1 is the Jacobi matrix between the global and
natural coordinate system. Given that the boundaries of the
element in the parametric space are determined by the
parameters
and
, mapping from the parametric to the
normalized space is also necessary.
1
1
det 22
j
j
ii
k
J (40)
The piezoelectric coupling matrix describes the
coupling between the mechanical and the electric field, and
it is given as:
TT
Td
uu u
V
V
KK BeB
(41)
7 NUMERICAL EXAMPLE
The considered two examples use a model of a
composite cylinder arc with two piezolayers. The structure
is a simply supported composite cylindrical arch of the
radius R = 100 mm and the width b = 62.8 mm. The
composite layers have a "balanced" stacking sequence [45/
− 45/0] s (orthotropic material properties) and the
piezolayers, which cover completely the top and bottom
surfaces, are oppositely polarized. The thickness of each
composite layer is 0.12mm and that of the piezolayers is
0.24 mm. The properties of the composite and piezoelectric
layers are given in Tab. 1 [33].
Deformation of the arch is achieved in two different
ways with purely mechanical loads and with loads obtained
through excitation of the piezolayers by electric voltage.
Although the symmetry of the structure allows to
consider its quarter with appropriate boundary conditions
applied, it was decided to use the full model as the structure
is relatively simple and thus the resulting FE model would
not be numerically demanding. The structure is discretized
using a NURBS patch with the basis function of degree 4
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
Tehnički vjesnik 30, 2(2023), 416-425 423
in arch direction and degree 2 along the cylinder axial
direction (width).
Table 1 Material properties
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
The mesh has 20 elements and 25 control points in arch
direction. Along the cylinder axial direction (width), the
number of elements is 4. The obtained results are compared
with the results obtained using the classical finite element
model discretized with the Reissnerr-Mindlin shell type
developed for modeling composite laminates with
piezolayers [33]. This is a full biquadratic shell element
with 9 nodes. In the case of classical element, three
different finite element meshes were used to check the
convergence the selected number of elements in the
circumferential direction was 10, 20 and 40, respectively.
To check the suitability of the developed FE
formulation for both purely mechanical and coupled
piezoelectric cases, two different loading cases were
considered. The first one is purely mechanical. In this case,
the structure is loaded with a vertical force (1N/m)
uniformly distributed over the width line that splits the arch
into halves. In the second case, the excitation is achieved
by means of electric voltage of 100 V supplied to the
piezolayers. As the piezolayers have symmetric position
with respect to the mid-surface of the arch and are
oppositely polarized, this excitation induces bending
moments uniformly distributed over the supported edges
and the free cylindrical edges of the arch. The presented
results of the isogeometric approach indicate excellent
agreement with the results reported in previous study
(Fig. 6 and Fig. 7), but it should be noted that the
isogeometric formulation uses less degrees of freedom to
achieve the same accuracy as the classical FE formulation.
Figure 6 Simply supported composite cylindrical arch under uniformly distributed vertical force
Figure 7 Simply supported composite cylindrical arch subjected to excitation of the piezolayers by electric voltage
8 CONCLUSIONS
Active thin-walled structures have received a great
deal of attention over the previous two decades, which is
due to the obvious advantages they offer over the classical,
passive structures, recognized in the improved dynamical
behavior, safety and robustness. Suitable numerical tools
that enable their efficient modelling and simulation of their
behavior are an important prerequisite for their successful
and efficient development.
This paper proposes an isogeometric FE formulation
for shell structures made of composite laminates with
embedded piezoelectric layers. The piezoelectric layers
characterized by the coupled electro-mechanical field may
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
424 Technical Gazette 30, 2(2023), 416-425
be used both as actuators and sensors and thus enable active
behavior of the considered structures. By coupling sensors
to actuators via a controller, adaptive structural behavior is
obtained. The developed formulation offers a tool for
modeling mechanical and electrical field as well as their
coupling in such a structure, and it could be extended in a
further step to include a control algorithm and thus test its
effects. The formulation covers laminated structures. In
order to keep the level of numerical effort in an acceptable
realm, the formulation relies on a first-order theory with
the Reissner-Mindlin kinematics. Hence, the transverse
shear effects are included, which is important for laminated
structures made of composite materials, due to their
susceptibility to those effects.
The problem of defining a normal to the reference
surface was given particular attention. While it is a
relatively trivial task in the framework of the classical FE
formulations, it turns out to be a relatively cumbersome
task in an isogeometric formulation. In this work, the
approach based on the Greville's points is adopted.
The structure considered in the example and subjected
to both purely mechanical and electrical loading, shows
that the developed isogeometric formulation with higher
degrees of basic functions offers a good match for the
classical FE formulation of a piezolaminated shell. The
results are easily matched with less degrees of freedom and
therewith with a smaller numerical effort. This is quite
promising for further work, which will aim at extension of
the formulation for geometrically nonlinear and dynamic
analysis.
Acknowledgements
This research was financially supported by the
Ministry of Education, Science and Technological
Development of the Republic of Serbia (Contract No. 451-
03-9/2021-14/200109) and by the Science Fund of the
Republic of Serbia (Serbian Science and Diaspora
Collaboration Program, Grant No. 6497585).
9 REFERENCES
[1] Hashemi, K. S. H. (2021). Nonlinear vibration response of
piezoelectric nanosensor: influences of surface/interface
effects. Facta Universitatis-Series Mechanical Engineering.
[2] 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), 1130-1137.
https://doi.org/10.22055/jacm.2020.31899.1987
[3] 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),
61-76. https://doi.org/10.31181/rme200101061t
[4] Gabbert, U. & Tzou, H. S. (2000). Smart Structures and
Structronic Systems. Kluwer Academic Publishers,
Dordrecht, Boston, Amsterdam.
[5] 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. https://doi.org/10.1088/1361-665X/aad480
[6] 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), 372-381.
[7] Noll, M., Lentz, L., & von Wagner, U. (2019). On the
discretization of a bistable cantilever beam with application
to energy harvesting. Facta Universitatis-Series Mechanical
Engineering, 17(2), 125-139.
https://doi.org/10.22190/FUME190301031N
[8] Mohammadkhah, M., Marinkovic, D., Zehn, M., & Checa,
S. (2019). A review on computer modeling of bone
piezoelectricity and its application to bone adaptation and
regeneration. Bone, 127, 544-555.
https://doi.org/10.1016/j.bone.2019.07.024
[9] 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),
109-126. https://doi.org/10.22190/FUME210311036N
[10] 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, 11-
22. https://doi.org/10.1016/j.finel.2011.11.006
[11] Reddy, J. N. (2003). Mechanics of Laminated Composite
Plates and Shells, Theory and Analysis. Second Edition.
CRC Press. https://doi.org/10.1201/b12409
[12] Gabbert, U., Koppe, H., Seeger, F., et al. (2002). Modelling
of smart composite shell structures. Journal of Theoretical
and Applied Mechanics, 3(40), 575-593.
[13] Marinkovic, D. & Rama, G. (2017). Co-rotational shell
element for numerical analysis of laminated piezoelectric
composite structures. Composites Part B: Engineering, 125,
144-156. https://doi.org/10.1016/j.compositesb.2017.05.061
[14] Rama, G., Marinkovic, D. Z., & Zehn, M. W. (2017). Linear
shell elements for active piezoelectric laminates. Smart
Structures and Systems, 20(6), 729-737.
[15] Kulkarni, S. A. & Bajoria, K. M. (2003). Finite element
modeling of smart plates/shells using higher-order shear
deformation theory. Composite Structures, 62(1), 41-50.
https://doi.org/10.1016/S0263-8223(03)00082-5
[16] Balamurugan, V. & Narayanan, S. (2009). Multilayer
higher-order piezolaminated smart composite shell finite
element and its application to active vibration control.
Journal of Intelligent Material Systems and Structures,
20(4), 425-441. https://doi.org/10.1177/1045389X08095269
[17] 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), 269-
283. https://doi.org/10.22190/FUME190530030M
[18] 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), 345-357.
https://doi.org/10.1177/1045389X17705538
[19] 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), 4135-
4195. https://doi.org/10.1016/j.cma.2004.10.008
[20] Kiendl, J., Bletzinger, K. U., Linhard, J., & Wüchner, R.
(2009). Isogeometric shell analysis with Kirchhoff-Love
elements. Computer Methods in Applied Mechanics and
Engineering, 198(49-52), 3902-3914.
https://doi.org/10.1016/j.cma.2009.08.013
[21] Nguyen-Thanha, N., Valizadehb, N., Nguyenc, M. N.,
Nguyen-Xuand, H., Zhuange, X., Areiasf, P., Zih, G.,
Bazilevsg, Y., Lorenzisa, De L., & Rabczukb, T. (2015). An
extended isogeometric thin shell analysis based on
Predrag MILIĆ et al.: Reissner-Mindlin Based Isogeometric Finite Element Formulation for Piezoelectric Active Laminated Shells
Tehnički vjesnik 30, 2(2023), 416-425 425
Kirchhoff-Love theory. Computer Methods in Applied
Mechanics and Engineering, 284, 265-291.
https://doi.org/10.1016/j.cma.2014.08.025
[22] Milić, P. & Marinković, D. (2015). Isogeometric FE analysis
of complex thin-walled structures. Transactions of
FAMENA, 39(1), 15-26.
[23] 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), 276-289.
https://doi.org/10.1016/j.cma.2009.05.011
[24] Echter, R., Oesterle, B., & Bischoff, M. (2013). A hierarchic
family of isogeometric shell finite elements. Computer
Methods in Applied Mechanics and Engineering, 254, 170-
180. https://doi.org/10.1016/j.cma.2012.10.018
[25] Yujie, G. & Martin, R. (2015). A layerwise isogeometric
approach for NURBS-derived laminate composite shells.
Composite Structures, 124, 300-309.
https://doi.org/10.1016/j.compstruct.2015.01.012
[26] 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), 159-179.
https://doi.org/10.1002/nme.4730
[27] Dornisch, W., Klinkel, S., & Simeon, B. (2013).
IsogeometricReissner-Mindlin shell analysis with exactly
calculated director vectors. Computer Methods in Applied
Mechanics and Engineering, 253, 491-504.
https://doi.org/10.1016/j.cma.2012.09.010
[28] 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, 106-137.
https://doi.org/10.1016/j.cma.2014.07.020
[29] Piegl, L. & Tiller, W. (1995). The NURBS Book. Springer.
https://doi.org/10.1007/978-3-642-97385-7
[30] Rogers, D. F. (2001). An Introduction to NURBS with
Historical Perspective. Academic Press.
https://doi.org/10.1016/B978-1-55860-669-2.X5000-3
[31] Gerald, F. (1993). Curves and surfaces for computer aided
geometric design: a practical guide. 3rd. edition. Academic
Press Inc.
[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), 503-513.
https://doi.org/10.1177/1045389X06067139
[33] Marinković, D., Köppe, H., & Gabbert, U. (2006).
Numerically efficient finite element formulation for
modeling active composite laminates. Mechanics of
Advanced Materials and Structures, 13(5), 379-392.
https://doi.org/10.1080/15376490600777624
[34] Zienkiewicz, O. C. & Taylor, R. L. (1989). The finite element
method - Volume 1: Basic formulation and linear problems.
McGraw-Hill, New York.
Contact information:
Predrag MILIĆ, PhD, Assistant Professor
University of Niš,
Faculty of Mechanical Engineering in Niš,
Aleksandra Medvedeva 14, 18000 Niš, Serbia
E-mail: [email protected]sl
Dragan MARINKOVIĆ, PhD, Full Professor
(Corresponding author)
University of Niš,
Faculty of Mechanical Engineering in Niš,
Aleksandra Medvedeva 14, 18000 Niš, Serbia
Lecturer and Research Assistant,
Department of Structural Analysis, Berlin Institute of Technology,
Strasse des 17, Juni 135, 10623 Berlin, Germany
E-mail: [email protected]
Sandra KLINGE, PhD, Full Professor
Department of Structural Analysis,
Berlin Institute of Technology,
Strasse des 17, Juni 135, 10623 Berlin, Germany
E-mail: [email protected]
Žarko ĆOJBAŠIĆ, PhD, Full Professor
University of Niš,
Faculty of Mechanical Engineering in Niš,
Aleksandra Medvedeva 14, 18000 Niš, Serbia
E-mail: [email protected]
