Smart Structures

Professor Mohammed Rabius Sunny

Department of Aerospace Engineering
Indian Institute of Technology, Kharagpur
Week 07
Lecture No: 40
Analysis of composite laminate with piezoelectric patches – computer programming
Part 05

In this lecture, we will discuss about solving the actuation problem in composites with the
help of Piezos using MATLAB. So, we will write computer codes and we will see how this
analysis can be done. We have seen that this kind of analysis involves a lot of computation.
So, this is not something which is to be done manually for even slightly complicated case.
So, for that we need to develop a computer code and solve the problem and that is what,
we are going to discuss today.

So, this is the example structure that is going that we are going to take.

It is a composite with a set of layers, and it has two piezoelectric patches and at the top and
bottom. So, the top piezoelectric patch and the bottom piezoelectric patch are throughout
the surface of it. And in this example, there are six layers of composites 1, 2, 3, 4, 5, 6 and
at the two sides of it, there are two piezoelectric patches. So, there are total eight layers.
Among those eight layers, six are composite plies, and at the at the two opposite surfaces
these are piezoelectric patches.

These are axis: x axis, this is y axis. From the x axis, let us denote the dimension as a and
along the y axis, let us denote the dimension as b. And then, the piezoelectric patches would
be excited with voltage and we will see the response.

(Refer Slide Time: 1:58)

So, for that we have license from MATLAB for doing this exercise and the license is
available in the course portal here. So, there is MATLAB menu and under that there is a
MATLAB access and MATLAB training.

(Refer Slide Time: 2:30)

So, those who do not have much experience above with MATLAB they can check this
training and however, we are not going use very advanced features of MATLAB. Basic
features would do. Now to access MATLAB we can go to MATLAB here and the first two
links through those, the sign up can be done, and after signing up we can access the
MATLAB. So, after signing up we can come to this use link and that will take us to
MATLAB.

So, if I go to this use link, then it will give us two options: one is open MATLAB online,
one is installing MATLAB by going to the installation option, we can install it and use it.
And by using this option open MATLAB online, we can do it through the browser. So,
here I am doing it through the browser. So, if I open it here, it shows me this. So, the entire
MATLAB front end is in our browser now. And these are the course which has been
written. Now I am going to explain that.

(Refer Slide Time: 3:34)

So, here we are going to use a symbolic toolbox. So, we will use the symbolic features.
Initially we should do clear all, so in case there is something in any variable predefined
that gets deleted otherwise, that can create some problem. And then x and y are our
symbolic variables, which we are defining here. Now let us define the set of inputs.

So, here we have to first input the geometrical features. So, as we have seen that the plate
has length of a and width of b, length means dimension along x axis and width means
dimension along y axis. So, a is 0.5 and b is 0.3, now there are total 8 layers including the
pieces and composites. So, that we are defining here as a variable NN. So, it is NN equal
to 8.

Now, here we are writing the code in a more generalized way. So, we are considering the
fact that the top piezo and the bottom piezo can be may not be identical, in terms of their
thickness or in terms of their material properties. So, we define the variables Ectop that
means, elastic modulus of the top piezo.

So, anything with a subscript top means that it denotes the top piezoelectric patch. So,
Ectop means elastic modulus of the top piezoelectric patch, neutop means the Poisson's
ratio of the top piezoelectric patch. rhoctop means the density of the top piezoelectric patch,
and d31top means the d31 constant of the top piezoelectric patch. Similarly, we have
Ecbottom, neubottom rhocbottom and d31bottom, they are the same quantities for the
bottom piezoelectric patches. Now here, the quantities that we have given are same, but we
can keep it different and that is fine.

Now, comes the properties of the composite plies. Now in here, we consider the fact that
each ply has same material property and same thickness, but they can differ by – in fact,
the thickness not need not be same, but we are considering same material property and their
difference is in the orientation of the fibers. So, E1 E2 are the elastic modulus of each
composite ply along two directions. E1 is in one direction which means the fiber direction
and E2 is in the perpendicular direction. These are all unidirectional plies. neu12 is the
major Poisson's ratio and then if I apply the formula, neu21 is equal to E2 multiplied by
neu12 by E1, that gives me the minor Poisson's ratio. G12 is the shear modulus; rho is the
composite ply density and theta means the orientation of fiber in each of the plies. So, there
are 8 plies. So, the it goes on like this: 0, 90, 0, 0, 90, 0.

Again, if we go to the figure there is the first piezoelectric layer. So, here the orientation
is 0, then 90 then 0 then again 0, 90, 0. So, it is a it is an example of cross ply which is
symmetric, but again that is not a requirement, any kind of lamination sequence that we
give symmetric anti symmetric or whatever, it would work even it did not be cross ply, it
can be angle ply also.

Now, let us define the z's. So, z here means the z coordinate of each of these two sides of
the layers. So, at the bottom most surface we have z1. So, here z1 means distance of the
bottom most surface from the x y plane with the sign taken care of. So, bottom surface
means we are going along the negative z direction. So, that is a negative quantity. The
highest negative quantity in magnitude and then the next junction that is z2, then the next
one z3, z4, z5, it goes all the way up to z9 because, there are total 8 layers. So, there are
total 9 z quantities. So, if I subtract z1 from z2 that gives me the thickness of the bottom
piezoelectric patch. Similarly, if I subtract zNN plus 1 that means, the last z to the 1, but
last z. So, zNN plus 1, minus zNN, that gives me the thickness of the top piezoelectric
patch and then, if I subtract zNN, I mean, if I subtract z2 from zNN that gives me the total
thickness of the composite plies, that laminar thickness is – thickness of the entire
composite laminate is that.
So, these are z's are arranged. So, we go from bottom to top. Now, let us give the force
input. So, here the force can be distributed forces along the x y or z direction over the
surface or there can be voltage actuation. So, qx means if there is any distributed force
along the x direction that is q. So, here we are assuming it to be 0, we can give a constant
value here or we can write it as a function of x and y also. Similarly, qy and qz vtop means
the voltage that is applied to the top and vbottom means voltage that is applied at the
bottom.

Now, after analysis we might be interested in finding out stress, strain, or any such quantity
in any particular point and that is defined here. xp yp means the coordinate of that particular
point and Lp means at which layer. Now, here the layer means we are not considering the
piezo layer, here the layer is which composite layer. Again, if we go to the figure. So, xp
yp denotes the location in the x y plane.

So, x means the value of the x coordinate and y is the y coordinate and then comes Lp, Lp
means layer. So, while defining Lp, we call as the first layer and this we call as the second
layer. So, that count of the Lp layer starts from the first ply layer not from the – So, the
count for, Lp layer starts from the first ply, not the first piezo layer. So, it is 1 2 3 4, it is
goes up like that.

So, in our case Lp means, 2 means, it is this layer. First there is a piezo layer, then ply
number 1, then ply number 2. After that we have to go to define the Q matrix. Now, Q for
the piezo layers is defined considering them to be transversely isotropic. So, in the x y
plane, they are isotropic. So, if I know the Ec, I just divide Ec by 1 minus nu square, and
then I multiply the matrix 1 nu 0; nu 1 0; 0 0 1 minus nu by 2, that gives me the Q for the
piezo layer. So, Q 1 is for the first piezo layer and QNN is for the last piezolayer.

So, here we can see that we have defined Q as a 3 dimensional vector, the reason is that
the first 2 indices are for the Q itself, I mean it is Q is a matrix. So, the element goes here
that is why, we need the first 2 index and the last one, the last index is for denoting which
layer it is. So, this denotes the Q matrix for the first layer. Similarly, this denotes the Q
matrix for the last layer.

Now, we have to find out the Q matrix for the plies.

Now, while finding out the Q matrix for the plies because we are starting from E1 E2 nu12
nu21 and the G. So, first we need to find out the Q matrix with respect to the material
coordinate and then we have to transform it to the x y coordinate. Now, because the
properties of each ply are same. So, their Q matrix with respect to the 1 2 system is same.
So, that is why we just define it once and we get it. This is the method for that. This is the
equation for that, that we have discussed. And then we go through a loop which goes from
2 to NN minus 1, because the layer number 2 means the first ply layer, layer number NN
minus 1 means the last ply layer. So, that is how we traverse. And while doing that, we
define the transformation matrix. Here theta is the variable. So, theta denotes the
orientation of the fiber. Here it is i minus 1 because i starts from 2. So, when i is equal to
2, I am in the first ply layer. So, that is why i minus 1, it denotes 1.

So, then we apply the same formula that we discussed T inverse Q into T and that gives us
Qxy. And then, we store that Qxy in the 3 dimensional vector for Q. And accordingly, we
finally get a 3 dimensional vector that contains the Q matrices for each layer. Now, that we
have got the Q matrices, we have to find out the A B and D matrix. We know, the relations
for A B D matrix and that is what is done here. And again, it has to be added over all the
loops. So, we initialize A B and D to 0 and then, we write the formula for the A B and D
matrix. And finally, at the end of the loop the A B D matrix is fully calculated.

Now, we have to define the basis functions, that we are going to use for this analysis. So,
let us go back to the figure once again.

Here, the consideration that we are doing in this example, it can be changed also that is up
to the user. Here, the consideration is that all the 4 edges are simply supported and that tells
us that at x equal to 0 and a, and y is equal to 0 and b, which means at all the 4 edges my
w is 0. And the bending moment may be 0, bending moment should be 0. So, at x equal to
0 and a, Mx is 0 and that tells us that Dxx is equal to d 2 w by dx 2, plus Dxy plus d 2 w by
dy 2, is 0. And similarly, at y equal to 0 and y equal to b, My equal to Dxy equal to d 2 w
by dx 2 plus Dyy multiplied by del 2 w by del y 2 equal to 0.

So, our approximation should satisfy at least this boundary conditions. If I can satisfy this,
that is good, otherwise, even satisfaction of this, this would also serve the purpose. And
then we have at x equal to 0 and a, we have, V0 equal to 0. And we have Nx which means
Axx multiplied by del u0 by del x, plus Axy which means del V0 by del y, equal to 0 and at
y equal to 0 and y equal to b, we have u0 equal to 0. And the natural boundary condition is
Ny equal to Axy multiplied by del u0 by del x, plus Ayy del v0 by del y, equal to 0. So, these
are the set of essential and natural boundary conditions, if we can satisfy this in natural
boundary condition also that is good. So, in this problem the function that we have chosen
satisfy all the boundary conditions.

(Refer Slide Time: 17:15)

So, this these are what we have as phiu, phiv and phiw.

So, phi u is sin pi x by a, cosine pi y by b, and sin 2 pi x by a, cosine 2 pi y by b. Now, here

we have just chosen 2 terms, we could have chosen multiple terms. In fact, while doing a
problem in its entirety, we should keep increasing the number of terms until, the
convergence is achieved. Here, for the sake of simplicity just for the demonstration purpose
we are showing it with 2 terms. We could have chosen terms like sin pi x by a, multiplied
by cosine 2 pi y by b, and also sin 2 pi x by a, multiplied by cosine pi y b, and so on.

Then there is a variable Nu, Nu is the length of this vector phiu. So, this basis functions are
put in this array phiu and that is the length of it and that is Nu. Now, while discussing the
theory we denoted this variable as m, but here this m might be used in for mass matrix also.
So, that is why to avoid any conflict of variable names, we are using Nu here. Similarly,
phiu is here, phiv is here. So, we are using this functions to define phiv and Nv is the length
of it. And phiw, these are the functions used to define phi w and these are the length of it.
So, it can be easily verified that this basis function satisfies all the boundary conditions.

Now, here we are defining the derivatives, because we are using the symbolic calculation.
So, we can easily do the derivatives differentiation symbolically. And in MATLAB, the
function that is used to do the differentiation is diff, d i f f.

So, diff phi u means which variable I am trying to differentiate and x means with respect
to which variable. So, phiux means derivative of phiu with respect to x. So, phiu has two
constituents, phiu1 and phiu2, this is my phiu1 and this is phiu2. So, it differentiates both
and puts in an array phiux. Then we are finding out the derivative of phiu with respect to
y. And then here we are differentiating phiv with respect to x. And then here we are
differentiating phiv with respect to y. Now, here we are differentiating phiw with respect
to x, and putting it in a variable in an array phiwx. Here we are differentiating phiw with
respect to y.

Next, we need to calculate the second order derivatives of phiw. So, that is what we are
doing here. Again, under differentiation phiw is what we are differentiating x is with
respect to what and 2 is the order of differentiation. So, when we find out the derivative of
first order, we did not need to put the order, but if we are doing it for second order or any
higher order, we need to put the order and that is what we are doing here. So, this is
derivative of w with respect to x done twice. This is second order derivative of w with
respect to y. And it is a mixed derivative, derivative of w with respect to x and y. So, this
is one derivative. And again, this entire thing is put under another derivative. So, that is
how the mixed derivative is done.

So, after having defined these derivatives of the basic functions, now we define the
quantities like NI, NW, BI, BW because if we go back to the formulation, we saw that our
final equation contains these matrices like MII, MIW, KII, KIW and so on.

(Refer Slide Time: 21:38)

And to find out these, we need these quantities NI, mII, NW, BI, BW and so on. And that is
what we are doing here.

(Refer Slide Time: 21:52)

So, this is NI. So, as per the definition that we discussed while discussing the theory we
have put it. So, in the definition of NI, first is phiu then 0s. Then in the next we put 0s, then
phiv. And this is the definition of NW that we discussed. Now, we define BI which involves
derivatives of the basis functions, and that is how it is defined. So, again it is a row and
below it is another row and below it is another row and it is Bw. It has 3 rows, the
derivatives of phi is of different order with respect to different variables.

Now, we need to find out the thickness of the composite layers. So, total thickness of all
the composite layers. So, we have seen that – if we again go back to the image, the total
thickness of the composite layer is basically the difference of z at this point, minus at this
point. So, at this point, the z is zNN minus 1, and here z is z2. So, if I subtract z2 from zNN
minus 1, that gives me the total thickness of all the composite layers, and that is what is
done here. And then we have found out the thickness of that bottom piezoelectric patch
which is z2 minus z1, and the top piezoelectric patch which is zNN plus 1, minus zNN.

Next, we find out the quantities that are related to inertia. So, we have defined m, S and I,
which are the mass, the first moment of inertia, and second moment of inertia and
accordingly they are put in mII, mIW, mWI and mWW matrix.

Next, we find the components of the mass matrix. So, again if we go to the final equation
that we use – here we have seen that, there are 4 matrices M II, MIW, MWI and MWW, and the
formulas are like this. So, applying those formula, we found out the MII, MIW MWI and
MWW matrix. Now, here one thing to note here, we are doing it symbolically. So, we can
do symbolic integration and put the limits using the MATLAB symbolic toolbox, but that
is not a general practice. In general, these integrations are done using numerical techniques.
For example, while using the finite element method, these integrations are done generally
using the Gauss quadrature formula.

So, here symbolically the integrations are done. Now these are double integrals because
the domain is 2 dimensional. So, we do that integral along x from 0 to a, once and then
again, we evaluate the integral from y equal to 0 to b. And finally, this has been converted
to a double variable because it is being done symbolically it might show it as a fraction.
So, a double comment converts this entire result as a result to a double variable. So, we
find out the constituents of the mass matrix.

Here using the same process, we find out the constituents of the stiffness matrix. We
already know the formulas; we can just apply it and do it.

Now, comes Np and Mp the block force and block moments. Now these are the relations
that we used for block force and block moments. Now this relation for the block force was
for one piezoelectric patch. If I want to do it for 2 piezoelectric patches, we can do it
separately and add it up and that will give us the total block force for all the piezoelectric
patches. So, here it was considering one patch and here it again it was considering 2
patches, 2 identical patches. So, this considers 2 identical patches. If the patches are not
identical, then the formula would be little changed. So, it would be tc, suppose we are doing
it for tc bottom. So, it would be tc bottom. And here I would have tb by 2, plus tc bottom by
2, and same thing here. And accordingly, if I am doing it for the top piezoelectric patch.
So, it is tc top. Here it is tb by 2 plus tc top by 2. And then, we will do it for 2 patches and
add them up that gives us the total block moment.

(Refer Slide Time: 27:08)

So, here it has been done in the generic way considering the piezos can be separate, their
voltage can be separate, their properties can be separate. So, you find out the Np and Mp.
So, we do it for 2 patches and add it up and the added quantity is this. And similarly, the
added block moment is this. Having found out Np and Mp, we can find out the force vector.
Force vector has 2 parts, FI is the in plane force vector and FW is the out of plane force
vector and we are already familiar with the formula for the in plane and out of force vectors.
Considering, the block force, block moment and the distributed in plane and out of plane
forces and those also involves evaluating integrals in 2 dimensions and that is what we
have done here. Just we have applied the same integral formula symbolically.

And finally, after calculating the M, K and F matrix, we are assembling everything here.
So, that gives us one particular set of equations. So, here for example, if I combine
everything and write it in this way MII, KIW, MIW. So, combining everything it gives me
full mass matrix, full stiffness matrix, and full force vector. So, then that gives me equation
of the form Mq double dot, plus Kq equal to F. By solving, which we can find the entire q
vector which has qu, qv and qw.

Now, here we are solving it as a static problem. So, the mass matrix we do not need. So,
here we are just solving Kq equal to F. And there is a command in MATLAB that is a
linsolve. So, that solves the system of equations Kq equal to F, and the returns the variable
here. So, q is our solution.

Now, once we get our q, that means, we know all the unknown coefficients associated
with each of these basis functions. After that, we can find out our quantities of interest like
strain components and the stress components, at our point of interest and we defined our
point of interest as a point with coordinate xp yp and composite layer number Lp. So,
epsilon0p that means, epsilon 0 evaluated at that point p that is independent of the z
coordinate. So, we find out the value of epsilon 0 at that point by this. We know that if I
multiply the BI matrix with q, that would give me the epsilon 0. Now here, if I just multiply,
it will give me the symbolic variable. And then, if I want to put the value of xp and yp at x
and y, we have to use this comment subs under this everything. So, this expression and in
the expression, x y should be substitute by xp and yp. And accordingly, kappa is defined.
And after defining kappa, we have to find out the zp that means, the z coordinate of the
point.

And we are finding out at the mid location of the Lp th layer. So, for that we can just use
the formula and that will give me the z coordinate of the mid location of our layer of
interest. Then, we can apply this equation, epsilon equal to epsilon 0 minus z into kappa,
and that is going to give us the epsilon vector, the strain vector. And if, I multiply the strain
vector with the corresponding q vector because if the Lpth layer. So, Lpth composite layer.
So, absolutely it is 1 plus Lpth layer.
So, if I multiply the corresponding q matrix with that epsilon vector, that is giving us the
sigma, the stress at that layer. So, we have found out the stresses at that particular point.
So, accordingly we can be interested in find out stress strains at multiple points. Then, we
have to define multiple points and we have to put this entire thing under a loop.

Now, we might be interested to look at the deflected shape for that what we can do is we
can divide the entire surface into a grid point. And let us divide that into a 20 by 20 grid.
So, Ndivx means the number of divisions along the x, Ndivy means number of divisions
along y.

So, here we define a variable xn and yn. So, that will contain the x and y coordinate of the
of our grid points. So, if we draw it it looks like this. So, here I have the grid and by grid
points I mean these points. So, we are finding out the x and y coordinate of all of these
points.

(Refer Slide Time: 32:35)

And that is what we are doing by that piece of the code. So, here we are finding out the
coordinate of the grid points. And then at those grid points, we find out u, u0 v0 and w.
Again, we use the substitute command because we have found out our q.

So, if I multiply the first any elements of q with the components of phiu that gives me u 0.
And here, it will give me as a symbol. So, on on that if I substitute x and y by this grid
points are the coordinates, numerical values that would give me the actual value of u0 at
that point, and same thing we do for v0 and w. So, un means the value of u0 at the grid
points, v means the value of the v0 at the grid points, and wn means the value of the w at
the grid points. It is a 2 dimensional vector because i denotes the grid point number along
the along the x direction, j denotes the grid point number along the y direction. So, finally,
after that we can do a surface plot and we can visualize the displacement.

So, now after looking at the code if we run it. So, it is running now. If I want to stop it in
between, we can just tap the button it would stop the code. So, now, it has run and it is
showing how w is varying along x and y. Now, we have not put the x and y level that we
can put easily. So, if I type x level and then it is X the unit in meter.

So, it puts the xlevel here. Similarly, ylevel also can be put. So, that is the y level. So, it is
y direction in meter and then z level is the displacement w, that we were we are visualizing
that is zlevel. So, it is w.

So, now, accordingly we can plot u0, v0 also. So, if I say un and if I just evaluate that
particular line itself, that would plot u0 for me. If I say v here, if I evaluate that that would
plot v0 for me.

And also, we were interested in finding out the stress which has already been calculated.
So, if i write sigmap because, that is the variable that we used for the stress at that pth point,
point p. So, if just type sigmap that is going to give me the sigma p as an array, the first
quantity is our normal stress along x direction, second quantity is normal stress along y
direction, this is the shear stress in the x y plane. So, this was about the code considering it
is a static problem.

If, we want to consider it is a dynamic problem then, the entire code almost remains same.
Then, the entire code almost remains same, only thing is that at the end here mass matrix
would be used.

And at the end, we have to do the time marching using the beta new mark method. So,
again if we just quickly go to the beta new mark method that we discussed. So, this was
the beta new mark method considering a system Mq double dot plus Cq dot plus Kq equal
to F. And this is the time marching scheme that we discussed. So, at first, we find out q.
So, here n plus 1, means n plus 1th time step in our code that is I. So, if we want to follow
our code then in the beta new mark part n plus 1, is i and n, is i minus 1.

So, at the i th time step which is n plus 1 here. So, q dot is this, entire quantity which is a
function of the qn double dot, qn dot and qn at the previous time step which is known to
us. Now here, we have not considered damping. So, these quantities associated with C
matrix would be 0. And after finding out qi double dot, we can find out qi and qi dot using
these two relations.

(Refer Slide Time: 37:36)

So, that is what has been implemented in the code here. So, we divided the entire time
domain here, the entire time domain T was divided into set of time steps.

So, here the deltat is 0.005, that we have chosen it can be different depending on the
problem. And then, suppose we are running the analysis till 0.8 second. So, we write an
array like this that contains values like 0, 0.005, 0.01, 0.015 and so on till 0.8. That means,
it contains the values of all the time instance and that is put in an array ti and Nts is the
number of time steps, which is in this case that is going to be 161 because our deltat is
0.005 and our ti is 0.8. And also, the voltages can be a function of time. Similarly, the qx,
qy, qz also can be function of time. Here, we have assumed that they are function of time
of the form omega t. So, the voltage can be V sin omega t, qx can be q sin omega t, and so
on. So, for that we have defined the corresponding omegas. So, omegax means the
frequency corresponding to qx, omegay means frequency corresponding to qy, and omegaz
means frequency corresponding to q, sorry. So, here omegax means frequency
corresponding to px, omegay means frequency corresponding to py and omegaz means
frequency corresponding to pz and omegaqdd is the frequency corresponding to the
voltage. So, this voltage also can be a function of time. So, if it is sin omega t, the
corresponding omega is this.

So, what we do here is – we run a loop. So, initially we initialize the quantities q, qd and
qdd to 0, because at the initial time step, time t equal to 0, everything is 0. And, in beta
new mark method, we have some parameters beta and gamma, we have set the parameters
here. And then, we are running a loop from 2 to Nts number of time steps.

Now, here the force matrix has to be calculated at each and every time step, because that
is a function of time. So, we have multiplied the expression of the force matrix into a
quantity sin omega t, and the omegas whether it is omega x omega y and omega z that
depends on which component we are evaluating. So, here we are evaluating the it
corresponding to the voltage. So, it is omegaqdd and here it is corresponding to the force
along x and y. So, it is omega x and omega y.

Now, we after evaluating this FI, we have got FI and FW as a function of time. We get our
force vector at each and every time step. And then from that using the technique that we
described, we find out q double dot, qdd means second derivative of q with respect to time,
qd means first derivative of q with respect to time and q is q. So, evaluate that at each and
every time step and we do the time matching. And, each after evaluating at each time step,
we might be interested to find out u, v and w or quantities like epsilon and sigma, that we
discussed before at our points of interest at each and every time step. So, here we are
evaluating up, I mean, u0, v0 and w. So, here we are not evaluating u0, v0, w in the entire
grid, rather here we are evaluating only one point. We could have found out u0, v0, w in the
entire grid also as a function of time, but we are not doing that here to keep it simple. So,
we are evaluating u0, v0 and w at point p at every time step and that is what is done here
the same substitute command.

And here, we are finding out sigma and the time loop ends. So, if I run this code now. So,
here we are printing this i. So, if I look at the command prompt, I can know at which time
step the code is running. So, it is soon going to be at 160 first time step and come out of it.

Now, if I want to find out w as a function of time. So, we have seen that w is a one
dimensional and now because we are not evaluating w at each and every point, we are
evaluating at only one point. So, the only index that it needs is the time index I. So, if I
evaluate w as a function of time. So, ti is the time and corresponding w. Ok, it was wp not
w the variable. So, we can see w has been evaluated here. So, w is varying as a function
of time in this fashion.

Accordingly, we can plot the variation of sigma, variation of epsilon and variation of other
displacement quantities. And we can put the x y level also by using the command that we
saw or, by using these buttons xlevel, ylevel and so on.

So, that was an overview of how we can do these calculations by writing a simple code in
MATLAB. And again, the code has been written in a more generic way. So, the
piezoelectric patches can be different and the basis functions, that we have chosen, here
are sin functions depending on the boundary conditions, they can be different and
accordingly proper basis function has to be put.

So, with that I would like to conclude this demonstration here.

Thank you.