Tip:
Highlight text to annotate it
X
In
the previous lecture, we have talked about the classical plate theory.
As an example of a problem which leads to a fourth order partial differential equation
which we had said that if we work things out it will be this is in the area A where w as
the transverse deflection, D is the flexure rigidity constant and q is the applied distributed
load on the top or bottom surfaces of the body. This led to a variational form which
was in terms of the second derivatives of w. Since, it is in terms of variational form,
in terms of the second derivatives of w, we require these derivatives of w to be defined.
We need continuity of w del w del x and del w del y. So, this is what we had required
as far as continuity and we said this would lead to just analogous to what happened in
the Euler Bernoulli beam theory in the 1D case to C one continuity requirement in two
dimensions. We had said that we will look at certain ways of constructing these C one
continuous elements.
One way of doing it is let us say I have a rectangular domain, I will take only a few
elements; so, this is size a, this is size b, let this be our x and y direction, so these
are all rectangles. So what we want is at the interfaces of these rectangles with each
other including these end points, I should have continuity of del w del x, w and del
w del y. If I can construct basis functions which can satisfy continuity of these quantities
then we are in good shape. There were several attempts to make such elements and various
versions of triangular, here we are not considering triangular; triangular and rectangular elements
have been reported in the literature, various types of such elements. The major problem
here is to satisfy these constraints exactly in the general case.
Let us take the simple case and look at an example that let us say this is my element,
rectangular element. This is my rectangular element with the nodes x1 tau, y1 tau, here
x2 y2 and here is x4 tau and y4 tau . What is the basic idea in constructing such an
element? We say that w as a function of x and y is equal to, a0 plus a1x plus a2y plus
a3x square plus a4xy plus a5y square plus a6x cube plus a7x square y plus a8xy square
plus a9y cube plus a10x to the power of 4 plus a11x cube y plus a12x square y square
plus a13xy cube plus a14y to the power of 4. Why did I take this whole polynomial? Because
if I would like to have continuity of w, del w del x, del w del y at these 4 points of
these 4 nodes of the element. I will have to have w, del w del x, del w del y as the
variables at each of these nodes. We will call them w, this is essentially theta
y, theta y is the rotation about the y axis which is given by del w del x and theta x
is rotation about the x axis. I am not bothered, I am not taking here the positive or negative
signs that one can find out; this is equal to del w del x. This is del w del y and this
is del w del x. So these 3 quantities are specified at the each of these nodes. There
are 12 quantities which are specified at all the 4 nodes combined for this element. I need
to have a polynomial which has at least 12 coefficients. But, here we said if we look
up to the cubic in 2D, it has 10 coefficients - unknown coefficients. We have to do more
than that; so, we go to the fourth order. When we go to the full fourth order there
are total of 15 unknowns. We need only 12 coefficients, so what do we do? Out of the
fourth order we retain only the so called symmetric terms, this term and this term . While,
we drop off from the fourth order polynomial these terms. What we are left with is the
whole cubic plus now I will rename this, I will call it a10 here and I will call it as
a11. This is the polynomial that we are going to fit to these values.
But how do we do it? Again the idea is simple that our function w (x, y) in the element
is in terms of 12 coefficients such that the value of this function w at the nodes xi tau,
yi tau then thetay at the node xi tau, yi tau, thetax at the node xi tau, yi tau which
we are going to call by a new name; this I am going to call as wi, this is thetayi, this
is thetaxi. The value of this function equals the coefficients at these nodes. So we have
w at xi tau, yi tau is equal to a0 plus a1xi tau plus a2yi tau plus up to a11xi tau yi
tau cube. Similarly, thetayi at xi tau, yi tau is equal to del w del x at the node xi
tau yi tau. This would be equal to a1 plus, a2 would not be there, 2a3xi tau plus up to
a11yi tau cube.
Similarly, we can give thetaxi is equal to del w del y at the node xi tau, yi tau; so
this would be equal to a2 plus up to a11 into 3xi tau yi tau squared. We can write the thetaxi,
thetayi, and wi at each of these 4 nodes which have coordinates xi tau, yi tau in this expanded
form. How many values we have? We have 12 values; it is a twelfth order, there are terms
in this expression. We can write it as, this whole thing wi, w1 thetay1, thetax1, then
w2, thetay2, thetax2 up to all the way down to w4, thetay4, thetax4 this is equal to this
matrix A into the vector of this coefficients a0 to a11. What are the entries of this matrix
A? It is essentially one xi, yi depending on the expression for thetaxi, thetayi and
wi in terms of the xi's and yi's, so then implies we can get a is equal to A inverse
into y where this y is this vector. I can find the coefficients of this polynomial expression
for w in terms of the values of w, thetax, and thetay and the nodes and all we can do
is essentially, define the finite element solution in terms of these nodal values of
w thetaxi and thetayi. This expression can be easily obtained once it is obtained then
this is known; now we go ahead and do the same job for all the elements and we have
the representation of the solution in the element and also globally.
This way we can construct wFE. What will happen is wFE is nonconforming in this case.
What do we mean by nonconforming is that, it does not strictly satisfy the constraint
of continuity of value and derivative on each of this faces on each of this edges of the
element from this side and from this side . The two elements sharing this edge will
have on this edge it can be shown continuity of w and the tangential derivative of w. Not
the normal derivative, the normal derivative here from the two sides
is not continuous. For example, on this face on this edge del w del y will not be continuous
but del w del x will be and w will be. This is the problem because, we required a weak
solution to have a minimum continuity condition of all these derivatives and w but we are
violating it, because of this violation finite element solution suffers. In many cases we
get pretty decent solutions and this can be used. Can we improve it? Yes we can improve
it and where does the improvement come from?
This is the rectangular element, so here this essentially the x direction, this is the y
direction of the element. In this direction in the x direction, what if I took this hermite
cubic function as we had created in the case of the beam analysis? If I take this one this
is what we had as N1 hermite cubic as a function of x. This one would be N2 hermite cubic as
a function of x. Similarly, I can create the other hermite cubics; this will be N3 and
N4, this will be N3 hermite cubic as a function of x and this will be N4 hermite cubic as
a function of x. What does the hermite cubic do? If I have
another element here and I have the same hermite cubic here also same definition, it ensures
continuity of the derivative also in this direction which is essentially the x derivative
del w del x. Similarly, I will define now the hermite cubic with respect to y in this
direction. so the hermite cubic with respect to y would be this function. These functions,
so I will call this N1 hermite cubic with respect to y, this is N2 hermite cubic with
respect to y, this is N3 hermite cubic with respect to y and this N4 hermite cubic with
respect to y. These are our hermite cubics, so this is N4, N3, N2, and this is N1 so we have this hermite
cubic in the x and the y direction so in tensor product family, for the quadrilaterals same
thing we can do here. How many functions we can create now by taking products of these
hermite cubics in the x and y direction? We will create 16 functions because they are
4 in the x direction. There are 4 in the y direction and these functions if we remember,
will now be given in terms of nodal values of some quantities.
It turns out that it will be in terms of nodal values of w del w del x, del w del y and del
2w del x del y. These 4 nodal values in terms of these 4 nodal values I can define the hermite
cubic. One thing has happened; this also enforces an extra constraint that we are asking for
this cross derivative of w that del two w del x del y to be also continuous in this
case. Nevertheless our basic requirement is satisfied of this w del w del x, del w del
y being continuous; so how can I construct this hermite cubics? At every node we have
wi this we had said theta yi, this was thetaxi and this will call it has kappa x yi. Essentially,
we will get the values the representation of w in terms of these 4 values at the node
that is total of 16 values. It is very easy to construct this functions,
so I will call my N1 hermite cubic as a function of x and y will be equal to the first the
hermite cubic corresponding to this line and this node it is N1 hermite cubic due to x
and hermite cubic which is one here, value one corresponding to this line so it will
be N1 hermite cubic y. Similarly, N2 hermite cubic xy will be equal to, let us say I will
take the same value one here and the derivative one here. It will be N1 hermite cubic in this
direction x. Let me take derivative in this direction, value in this direction so I will
take N2 hermite cubic with x, N1 hermite cubic with y similarly, N3 hermite cubic with respect
to x y is equal to I will have N2 N1 hermite cubic with respect to x and N2 hermite cubic
with respect to y. N4 hermite cubic with respect to xy becomes N2 hermite cubic with respect
to x and N2 hermite cubic with respect to y. We have essentially the 2 hermite cubic
in each of these directions which are non zero which have some meaning here which are
non zero at this point. Either in the value or in the derivative with respect to those
hermite cubics, products of those hermite cubics we construct the 4 functions here.
N2x into N2y will correspond to this cross term that is product of the 2 derivatives.
While, the value here will correspond to N1h hermite cubic will correspond to the value
wi and the others will correspond to the slopes in one direction or the other direction.
I can say that this will correspond to wi, this will correspond to thetayi, this will
correspond to thetaxi, this will correspond to kappa x yi. At every node so I can define
this hermite cubics with respect to the non zero hermite cubics . I will have the N1,
N2, N3, N4 here then I will have N5, N6, N7, N8 in terms of the hermite cubic which are
active here and so on. This way I can create all the 16 hermite cubics and now I have so
called conforming approximation over the rectangular elements. Whether, for this kind of a domain
if I make a mesh of quadrilaterals, I make the mesh of the quadrilaterals something I
will do here. In this mesh of quadrilaterals whether, this hermite cubics do satisfy the
continuity of x and y derivative on these faces. If they do then they are actually conforming
in all cases. If they do not, then they are not but for the rectangular domain they will
satisfy conforming condition; we check here this is an assignment problem which you can
work out and check.
Plate theory: since we are in the topic of plate theory is that, we have dealt with now
the Kirchhoff plate theory was for very thin plates where the shear effects could be neglected.
Let us go little further, so we can have a plate theory which is called Reissner-Mindlin
plate theory which is an improvement of over the plate theory that we have discussed which
led to the fourth order differential equation. The idea behind this plate theory is this
kind of a representation of the displacement field in terms of the transverse coordinate.
Here I will have this has theta y, v (x, y) should be a function of z also, z is equal
to v0(x y) plus z theta x and w (x, y, z) remains w (x, y). What this theory does? It
assumes that ezz is still equal to 0. This means gamma xz, gamma yz are constant While,
from our standard strength of materials we know the gamma xz and gamma yz is actually
parabolic across the in the in terms of z. Nevertheless, this is an approximation which
is made. I am not going to dealt too much into this plate theory as such but let us
bring out an important aspect of it. What happened in the case of Kirchhoff plate
theory that we said these things were 0? Because these things were 0, this theta y and theta
xx came out as minus del w del x and del w del x and del w del y. These are essentially
free or independent variables so in this problem we solve for u0, where is u0, v0 coming from,
if I have in plane loading? So u0 v0, if I do not have in plane loading then we can ignore
these two parts where z is with respect to again the centre line of the plate. So u0,
v0, w, theta x and theta y these are the 5 unknown functions that have to be solved for
in order to construct the full solution. This is one feature that we have further Reissner-Mindlin
plate theory as compared to the Kirchhoff's theory where we only needed to solve for u0
v0 and w3 functions.
If I do this, then let us see now we should able to tell me that exx is actually equal
to u0 comma x plus z thetay,x, eyy is equal to v0,y plus z thetax,y, ezz is equal to 0
gammaxz is equal to w comma x plus thetay gammayz is equal to w,y plus thetax and gammaxy
is equal to u0 comma y plus v0,x plus z theta y comma y plus thetax, x. This is the state
of strain at any point given in terms of u0 v0 thetay and thetax and w. In this strain
terms the first derivatives of w u0 v0 and thetax thetay are sitting. Only first derivatives,
so if I now write the strain energy that is I write the stress in terms of the strain so
the stress will also are in terms of the first derivatives of all these quantities. The product
stress into the strain sigma xx into exx will be in terms of the first derivatives of u0
v0 and so on.
The strain energy definition as we had written strain energy is integral over the volume
of half of sigma ij e ij dv, ij going from 1 to 3. This expression would be product of
at most first derivatives of u0 v0 w thetax and thetay. It is quadratic in terms of the
first derivatives of all these unknown functions. If I look at the strain energy for the strain
energy to be finite we only want these to be defined. We only want this first derivative
of u0 v0 with respect to x and y thetax, thetay and w all with respect to x and y to be defined
which means that here if I want to now construct a basis function to do the approximation we
need to use only C zero. I can use the standard C zero elements that we had talked about in
a very detailed way earlier when we started the 2 dimension problems. Those elements those
basis functions like we had the linear quadratic, cubic, triangles the tensor product, quadrilaterals,
the serendipity quadrilaterals all those things can be used for the approximation of u0 v0
thetax thetay and w in terms of x and y. This problem can be solved very easily using the
machinery of the shape function generation and integration and all those things that
we had developed earlier. This way we can essentially use the tools available to us
in a judicious manner. This model it is not to say that this model will do better than
the Kirchhoff's plate model or it will do verse. This model has its inherent problems
which can be corrected to take care of things. This one should read the books and enough
of material is available.
There are also higher order plate models which go beyond this expression of u in terms of
z theta y z theta x and we add more terms in terms of higher powers of z. We would have
z square, we will have some other quantity which is a function of x and y plus here also
z square more. This way we can construct so called higher order shear deformable theories.
In those theories shear stress or strain will be non zero and the higher ones will also
give a non constant shear stress. We will use sufficient number of terms to get a parabolic
variation of shear which is when we take cubic terms that is celebrated higher order shear
deformable theory. What is known as, HSDT due to Reddy and lots of other people have
done this. Avery good reference is on this plate theory is papers by g n Reddy papers
by Tarun khan from IIT Bombay and from Szobo, Actis and Babushka on a different approach
to creating this plate models. I will not deal with those things in detail.
I think we are in a position to deal with all these plate models and the refinements,
if we need them and how to handle them? Let us go little ahead with what we are doing?
Let us now talk of a more interesting problem in common practice now a day, because computational
tools have improved tremendously. We can get dual processor machines with very
high speeds 3GHz -3.4 GHz sitting on a desktop, at a very affordable price even here.
In that case, the capabilities have improved and because of that our desire to do more
refines modeling of the physical problems. If we see the beams, the plates, the bars
these are all idealizations of a 3D situation, we do 3 dimensional analysis.
I will briefly touch upon 3 dimensional analysis and how it is done? The language the procedures
are direct analogous extension of what we had in the 1 and 2 dimensional cases. Let
us say that here we have a 3 dimensional domain, I will take a cubical domain for simplicity
let us say the volume of the domain is v and the surface is given as delta v. Over this
domain, I would like to find the response of this structure when let us say on this
phase, I am going to fix my displacement completely fix this phase and on this phase partially
I want to apply a transverse load. This looks like a bending problem; now if it is a long
slender member we will see that what we get out of the beam analysis will be closed. If
it is flat and both dimensions are similar plate analysis will do a good job but that
is for us to see. Here we would like to do an honest 3 dimension analysis with the load
is applied here. How do we solve this problem? First of all
we are talking of 3D elasticity as an example of a 3D problem and for 3D elasticity we would
like to develop the variational formulation and then suggest how to go and create the
basic functions or the shape functions which can be used to solve the problem. Some things
that we should keep in mind before we proceed, we see that here I deliberately applied the
load in only part of the domain. This we should keep in mind and the boundaries now of this
3 dimensional domain are surfaces that also we should remember, we will have surfaces,
surfaces will have edges and the corners. These are phases, these are edges and these
are corners.
Let us do to write the 3D equation of equilibrium for the theory of elasticity sigma xx,x plus
sigmaxy,y plus sigmaxz,z plus fx is equal 0. Similarly, I have sigmaxy,x plus sigmayy
, y plus sigmayz,z plus fy is equal to 0. Here I am assuming that the stress tensor
is symmetric that is sigmaxz is equal to sigmayz, sigmaxy is equal to sigmayx so on and sigmaxz,x
plus sigmayz,y plus sigma zz,z plus fz is equal to 0. This quantity I am going to call
as r1, this is r2, and this is r3 that is these are components of the residue vector.
As we have done in the 2 dimensional case, let us take w with components w1, w2, w3 as
an admissible displacement, this is an admissible displacement, virtual displacement. While,
we will have the vector u with components u1, u2, u3 when each of these components u1,
u2, u3, w1, w2, w3 are functions of x y z; this is the unknown displacement field.
What we had done earlier? We had said that we will take the weighted residual formulation
that is we will take the volume take r dotted with w and integrated over the volume because
r was zero this will come out to be zero. Against any w is any admissible virtual displacement,
admissible means it should satisfy the geometric constraints wherever u is specified. From
our, let me jump the gun and say what we want our w to do in the standard situation.
We have not yet done integration by parts; we will do it is the w2 should satisfy the
constraints this one and this one.
Out of this I am not going to write the long expression, I will do integration by parts
once. By doing integration by parts an end of getting this expression sigma xx into exx
due to w plus sigmayy into eyy due to w plus sigma zz into ezz due to w plus sigma xz gamma
xz due to w plus sigma yz gamma yz due to w plus sigma xy gamma xy due to w this whole
thing dv will be equal to here f1w1 plus f2w2 plus f3w3 plus integral over delta v, I will
explain what is this traction vector t dotted with w dA.
This is that traction t what is traction t given as, if we remember t1 is equal to sigma
xx nx plus sigma xy ny plus sigma xz nz t2 that is the component of a traction vector
in the x direction t2 is the component of the traction vector in the y direction this
will be nx sigma yy ny plus sigma yz nz and t3 is equal to sigma xz nx plus sigma yz ny
plus sigma zz nz this is standard, if we do this integration by parts and use the green
divergence gauss, green divergence theorem we will get this in terms of that expression.
We will get this, now we see what are these nx, ny, nz these are the components of the
unit outward normal on the area; so here lets come back to our problem.
The unit outward normal is the vector n here, in this case since these are all parallel
n will be equal to 0i plus 0j plus k so it has only the component nz which is equal to
one on this phase. But in the general align phase we will have all the 3 components nx,
ny and nz .we will have all the components of n. If we see that where does the t become
non zero only on this part of the boundary area and in this case we see that t1 is 0
because I have drawn these as vertical forces t2 is 0, here t is equal to vector t is equal
to minus t3 e3 unit vector in the third direction.
I will put this expression for the t in the expression here and I will essentially go
ahead from here now again coming back to what we drawn there.
Now I am done the weak formulation the one I had written in the last expression is the
weak formulation, the weak formulation tells us that sigma xx due to the u into exx due
the w so on, have to be defined in order for this integral on the left hand side to be
finite which means sigma xx is in terms of the first derivative of u and first derivative
of v and so on. Similarly, for ex, ey and so on is all in terms of the first derivative
of u, v and w with respect to x, y, z.
Our weak form now requires only the first derivatives of the vector u and vector v and
vector w to be defined which means that in this case, we need with respect to what we
have done as far as the approximation is concerned, we need to construct a C zero approximation
in three dimensions that is we only want the value of this functions u1, u2, u3 and similarly
for the test function w1, w2, w3 which is the virtual displacement to be continuous,
derivatives need only to be defined. We now define continuous shape functions first of
all, how do we discretize the domain? Let us take our domain of the previous figure
and let us make a mesh here.
How will I make the mesh? For the mesh I have to honor the boundaries of the load profile
boundaries of material differences, so wherever my material has a change that boundary will
have the nodes and the edges and the vertices nodes and phases of elements sitting there.
What kind of element should I need? In this case, life is a little simpler so I will make
this kind of a partition, very simple partition we see here I will make with the different
colour so that we can see the effect. Similarly, in the third direction I will make the partitions
like this. If we see, what we have done? We have constructed elements which are actually
cubical in shape okay elements which are cubical in shape. Well the sides could have different
lengths so these are called in the language of finite elements brick elements. There are
other types of elements we can make but I am not going to discuss them here.
This is the brick element, so let us do a very simple job of constructing a finite element
basis functions or here we will concentrate only on the element so the shape functions
which are C zero continuous now for the brick element. Let us say co ordinate system is
this is x, this is y, this is z, x y and z. I will call this for the element, let us say
this is my generic element, this is my generic element tau, this is the node 1 of the element
tau node 2, node 3, node 4, node 5, 6, 7 and 8 so this will have co ordinates x1 tau, y1
tau and z1 tau similarly all the other nodes so the mapping now we will take this and I
am going to go straight away to the master element this will map to the master element
which is the master cube. We make a master cube like this; this will
be my psi direction, this will be my eta direction and this will be my zeta direction such that
this point is at psi eta zeta minus 1. This point is at psi zeta eta 1 just like
we had the 1D element master element then the 2D master element; now we have the 3D
master element. So, this point will have essentially psi your zeta will be 1 psi will be plus1
zeta will be minus1 eta will be minus 1, zeta will be minus1 so this will be plus1. So here
my psi will be minus1 zeta will be psi will be minus1, zeta will be minus1, eta will be
minus1 and zeta will be plus1. So this way we can construct the 4 nodes. We want to define
the basic functions or the shape functions with respect to these master nodes such that
they are such that the shape function is defined with respect to this node is 1 here, 0 at
all other nodes, 0 at all other master nodes. How do we define? Very simple; we now take
tensor of product of the 1D shape functions that we had created C zero shape functions
in the direction. N1 here, N1 hat as a function of psi eta zeta will be equal to N2 with respect
to psi N1 double hat with respect to eta N1 double hat respect to zeta and that is all
so I will take the tense of product as if this is my psi direction minus minus1 to plus1,
this is the eta direction from minus minus1 to plus1 and this is my zeta direction from
minus1 to plus1 so my confusion was because my psi direction is positive in x direction.
Here this shape function was non zero with respect to psi with respect to eta, with respect
to eta this term will be non zero with respect to zeta, this term will be non zero, so for
this zeta, this is the first one this is N1 double hat eta this is N2 double hat psi this
is N1 double hat zeta. In terms of this 1D shape function we can
construct these 3D shape function and we see there will be 8 of them. They will satisfy
completeness linear independence and all those properties in this way. Putting them together
we can construct the global basis functions remember that the basic functions will be
piecing together these shape functions from all the elements which share this node and
then we can construct the 3 dimensional finite element solution for each of these components
u1, u2, and u3 solve using the weak formulation that we have done and we are now in a position
to solve the 3 dimension problem, but remember that here applying the boundary condition
is a tricky.
Because we have to do integral of the given loads over the area of the phases, in general
for the general 3 dimensional domains defining the domain meshing the surface of the domain,
meshing the interior of domain is not such an easy job. For that we need sophisticated
mesh generators which have to be used to construct the mesh on the surface of the outer surface
and the internal volume. With this I will stop my brief foray into the 3 dimension problem,
before which we have discussed the plate problem which was essentially in between the 2D problem
and the 3D problem. Before that we talked of the honest 2D problem which is the planar
stress problem and planar strain problem as well as the 1 dimensional heat conduction
problem. Next we are going to develop methods for a
different class of problems. Problems which relate to Eigen value problems at the continued
level. For example, the free vibration analysis or the buckling analysis, how do we solve
those problems using the finite element method?