 We have been discussing analysis of 2-dimensional continuum, so we will continue with that, so what we did in the previous lecture was we considered a triangular element, so this element has 3 nodes 1, 2, 3, and at each node there are 2 degrees of freedom, so this is a 3-noded element with 6 degrees of freedom, so the field variables are U and V which are functions of X, Y, and T, and we interpolate the field variable in terms of the nodal values, that means U is interpolated in terms of U1, U2, U3, and V is interpolated in terms of V1, V2, V3, and the interpolation functions are derived, so that for example the interpolation function N1 at this node will have value 1, and at these 2 nodes it will have value 0, and in this particular case we use linear interpolation functions. We also discussed rectangular plane stress element, again the field variables here are U and V, this element has 4 nodes 1, 2, 3, 4, and at each node there are 2 degrees of freedom, therefore it's a 8-noded, sorry 4-noded element with 8 degrees of freedom. Now to facilitate the development of the element we made a coordinate transformation, so if you recall the elements of stiffness and mass matrices for the element have to be evaluated by carrying out these integrations, this is stiffness and this is, this is stiffness matrix, this is a mass matrix, so to facilitate this integration we transform this rectangular region to a square region so that integration can be done from minus 1 to plus 1. Now in the examples that we have considered so far it has been possible for us to evaluate these integrals in closed form, but this may not be always possible, for example if you have a quadrilateral element or a triangular element with curved edges or a quadrilateral element with curved edges like this, evaluation of these integrals will not be straightforward, so it becomes necessary to deal with 2 complexities, one is a complexity associated with the geometry and the complexity associated with evaluation of these integrals, so these 2 issues we will try to take a look in today's class. So we will consider to start with a linear quadrilateral element, so this has 4 nodes this is a coordinate system U and V, X and Y and U and V are the displacement fields, U is along X and V is along Y and this is the quadrilateral element with nodes 1, 2, 3, 4 and the nodal coordinates are X1, Y1, X2, Y2, X3, Y3 and X4, Y4. Now we need to now approximate the field variable U and V, variables U and V within this domain. Now what we do is we introduce a transformation so that a point in XY here gets mapped to a new set of coordinates XI and ETA that we have to introduce now what these are and this coordinate transformation is implemented exclusively for the purpose of evaluating these integrals, we will see what transpires when we implement this transformation, the main idea is to introduce this transformation to facilitate evaluation of the integrals, the stiffness matrix is given by this and the mass matrix is given by this, now the transformation that we are doing is from XY is taken to psi ETA, therefore we can consider the transformation in the mathematical form as X is a function of XI and ETA, Y is a function of XI and ETA. Now following the way we represent the field variables, we write now for this X and Y this expression where what we are doing is X of psi, ETA is being evaluated in terms of its nodal values and interpolation functions, similarly Y is interpolated in a similar manner, so what we are going to do is we are going to use the same interpolation function that we use for representing the field variables to represent the geometry also, that is geometry of this transformation. Now if you recall the N1 into N3, N4 for this configuration we have shown that this is of this form and compactly for all the four indices 1, 2, 3, 4 we have obtained this in the previous lecture, so we will use this now here. Now let us consider the line 2, 3, and this line 2, 3, now we want to investigate what is the relationship between the two, so let us consider the line connecting X2, Y2 to X3, Y3, upon making the transformation X, X equal to psi, ETA and Y equal to Y of psi, ETA how does this line gets mapped, or alternatively we can ask if this is a transformation we are using, what happens to this line 2, 3 in this configuration, whichever way you look at we will find the answer to the required question. Now let us therefore consider the line XI equal to 1 connecting 1-1 and 1, 1, that means XI equal to 1 this line, now at one of the vertices the coordinates are 1, minus 1, so what happens to X there, X is X1 into N1 plus X2 into N2, X3, X4, now how does this shape functions behave, except this N2 all other shape functions will be 0, so this becomes X2, similarly you can write the value of Y XI, ETA that means where does this vertex 1, minus 1 in this coordinate system go in the original coordinate system, so X has gone to I mean this, okay let us see, X2 we have found out similarly Y will go to Y2, that means if I am considering 1-1 that means the second point here it is going to 2 here, similarly the other vertex 1, now 1 you can verify that this goes to X3, Y3, now between these two point does the variation remain linear or not, that is the next question we have to see, now consider the line XI equal to 1, which is that line, this line XI equal to 1, now X of psi, ETA is X1, N1, ETA and wherever XI is 1, XI is there I will write it as 1, so I get this expression, okay, I am simply putting for XI 1, then I will use the properties of these functions and I will be able to show that this is given by this expression, it's a linear function of ETA, similarly the Y coordinate also get mapped to another linear function of ETA and I get X and Y to be this, now I can eliminate ETA from this and get the equation in XY plane what it is, so if you do that I get Y- this equal to this- this, so what I am doing is this I am taking to the left side, this to the left side and dividing this 1 plus ETA gets cancelled this is what I get, now this is actually the equation of straight line passing through X2, Y2 and X3, Y3, so we conclude that line 2 3 in figure A is transformed to line 2 3 in figure B, that means 2 3 is mapped to 2 3 here, so we can summarize this, this is the line 2 3 which comes here, psi equal to 1, psi equal to 1 is this line that comes here, ETA equal to 1 is this line it comes here, psi equal to minus 1 is this line it comes here and 1 2 is ETA equal to minus 1 that comes here, so that means that this quadrilateral is getting mapped to this, through the transformation that we are proposing, so these transformations of coordinates that is X, psi, eta in terms of N and nodal values of X and nodal values of Y for Y, this is quite similar to the transformation on field variable, this is what I have been pointing out, this is a field variable, UJ of T are the nodal values which are the degrees of freedom in our finite element model, that is what we conclude from this is the same interpolation functions these NJs are being used to represent both the field variables as well as the geometric variables, the displacement field and the geometry of the element are getting represented through the same set of interpolation functions, the resulting element formulation is called isoparametric formulation, okay, fine. So what we have achieved is we have mapped a quadrilateral to a square, so when we are performing integration instead of performing over a quadrilateral which can be quite you know unwieldy we have now the need to situation of being able to integrate over a square region, but what about the integrand, how does that behave? So let us consider the evaluation of the mass and stiffness matrices by transforming X, Y to psi, eta coordinates, this transformation we will now implement, so we need to digress a bit now to understand what the problem is, so first is we will consider a function F of X, Y in Cartesian coordinates and using this transformation X, I, eta if we map X and Y to X, I and eta what is the rule of transformation, see this is the rule as you may be knowing but how does it really originate? To be able to do that let us consider a point, this point with a position vector R in XY plane, so that is R is I X plus J Y, now psi and X are functions, psi and eta are functions of X and Y, so now I want that is this is a psi and this is a eta some coordinates, okay, so now I want to deal with this situation, so what is dou R by dou psi? It is I dou X by dou psi plus J dou Y by dou Y by dou psi, where I and J, I, J, K are unit vectors along XY and Z respectively, similarly dou R by dou eta is I dou X by dou eta plus J dou Y by dou eta, AB, the line AB if you see is what is this, so this plus this is equal to this, this is a vector sum, so if you are looking at AB it is simply the difference between the two which is this, dou R by dou psi into D psi, similarly if you look AD it is the vector difference of this and this, position vector of D and position vector of A and that AD is this, now the area element DA is given by the triple product of these two position vectors and the K along this, this gives the area, this is a result from vector algebra so you have, you can recall that, and if we do that and now put substitute for dou R by dou psi and dou R by dou eta the quantity that we have just now derived and this itself can be expressed as a determinant, see we have shown just here dou R by dou psi and dou R by dou eta and this is what I am substituting here and implement this cross product, you should know the rule of cross product and the dot product, so if you do that you get this expression. Now this quantity inside this parenthesis can be written as a determinant and the matrix associated with this is the Jacobian matrix and this is the determinant of the Jacobian, so DA is therefore Jacobian of the transformation into D psi D eta, therefore I is integral F of XY DX DY A becomes now minus 1 to plus 1, now that is one first part of our digression, the second is the question of numerical integration, so you may be familiar with Trapezoidal rule, Simpson rule and Newton quotes formula etc., but we will consider what are known as Gauss quadrature rules and I will explain what the rational of those rules are, so let us consider for sake of discussion and evaluation of an integral I cap A to B of X DX, first let us transform X into a new variable U so that the limits of integration are from minus 1 to plus 1, so what I will do, I will substitute U as alpha X plus beta, alpha and beta are the two parameters which I don't know, so when we are at the lower limits I want this to be minus 1 half, so I want U to be minus 1 half when X is A and U to be plus 1 half when X is B, so that helps us to solve for alpha and beta and the transformation I am looking for is given by this, so if you substitute that into this equation I will get minus 1 half to plus half F of this U now B minus A DU, so this is equal to B minus A into I where I is this integral, so I call this F of B minus A U plus A plus B by 2 as phi of U DU, so we will focus our attention on discussing how to evaluate this integral, so now what I wish to do is I have this function phi of U and this is my U axis from minus half to plus half, I want to evaluate this function at some sets of points U1, U2, some UN and I want to attach a weight R1 to phi of U1, R2 to phi of U2 and be able to evaluate this integral, so I am looking for evaluation of I by using this representation, so the effort here will be spent in evaluating phi of U, so any integration method that we use if we should get the ideal integration scheme is get a good approximation for I with minimum number of evaluation of phi of U, or in other words if you evaluate say phi at say 10 points which is the best way that I can combine them, right, for any other combination you know you can use trapezoidal rule and things like that then you can evaluate R, but then you are constraining yourself to place UIs in a particular manner, so the weighting factor also can be optimized, now the question is therefore how to select N that means how many number of points, how to select these weights and where to place these points where I am going to evaluate these functions so that I get a good approximation to I, keeping this UI plus 1 minus UI as a constant that means you sample it at a constant step size that may not be the ideal situation, because that depends on how this function varies, okay, it may not lead to the best solution always, now assuming that phi of U is continuous between minus half and half, we write now a power series expansion, convergent power series expansion in this form, so in this form, now this minus half to plus of a phi of UDU, now for phi of UI substitute this summation, assume summation which is this, 0 to infinity AM, UM, DU and upon performing integral I get this, so I can evaluate these terms and in terms of the unknowns A0, A1, A2, etc. you can get an expansion like this for the integral, now, so this expansion we have got, now we have this power series representation for phi of U, now if I evaluate phi of U at UI this will be the representation, okay, now therefore now I will substitute this into the representation that we are using, so phi of UDU is approximately this, this is what is our proposition and for phi of UI I will write this expansion, this is what I get, so I get now a representation I will interchange the order of summation and I will get a summation I mean representation like this, so in longhand if you expand this R1 into this infinite number of terms, R2 into infinite number of terms so on and so forth, now this I know right hand side what it is and the left hand side I have, I am rearranging the terms now A0 I will collect, coefficient of A0 are collected that is R1 plus R2 plus RN, A1 is this, M is this and this continues, now I will now compare the coefficients of A0, A1, A2, A3 on both sides I will get a set of equations now, for this to be equal coefficient of A0 here is 1 and coefficient of A0 here is this, this should be satisfied and similarly coefficient of 8, this is A0, A1 there is 0 here, so that must be 0, so by writing this we can write as many equations as you want, suppose if you are, if you consider 2N equation there will be 2N unknowns, what are the 2N unknowns? R1, R2, RN first set then U1, U2, UN, so that constitute 2N unknowns, now you can, you select these 2N equations and solve for those 2N unknowns and that gives you the representation that you are looking for, but how do you solve for that? This is a set of nonlinear equations, you can see quadratic terms here, so they are nonlinear equations, but there is a nice result which shows that if phi of U is a polynomial of degree not higher than 2N minus 1 then this U1, U2, UN which are the points where I want to evaluate the function according to the prescription that I just outlined would turn out to be the zeros of Legendre polynomials PN of U, that means these U's will satisfy this equation, so the roots and these roots are real, so once these roots are determined if I know U1, U2, UN then they are remaining N unknowns which are R1, R2, RN can be solved by solving a linear problem, is it okay? So I can find this constant as I wish, so if you are dealing with a polynomial then if you retain adequate number of terms you will be evaluating the integral exactly, there is no approximation at all, okay. Now for example N equal to 3, we have, this is a rule according to which we want to evaluate, so if you expand this now for N equal to 3 and carry out this differentiation I get this equation, so it turns out that the roots of the equations are U equal to 0 and plus minus half square root 3 by 5, so the first point I have to select is this, second one is 0, third point is this, with these weights if I add and evaluate the integral, so 2N minus 1 means 4, 6 minus 1, the 5th order polynomial will be exactly evaluated by these weights, okay, because a function is polynomial. Now I have done a few calculations here for different values of N the roots, these are the roots of the legendary polynomial, this is for N equal to 1, 2, 3, 4 and these are the weights, okay, so if you are implementing this you can use this information to do this. Now we can see how this formula perform, now if I want to evaluate integral minus 1 to plus 1 the answer is 2, so if I take one term it gives the exact answer, so the rule is a polynomial of order P is integrated exactly by employing N number of terms which is smallest integer greater than 0.5 P plus 1, so now X so 1, so P plus 1 is 2 you need at least one term, if you use one term you will get exact solution, but this is of course 0 because it's an order function under symmetric limits there is no problem, but you come here P equal to 2, so 2 plus 1 is 3 1.5 I need at least 2 terms before I can get the value of the integral correctly, so if you see here if I take only one point I get answer as 0 which is wrong, if I get 2 points this is the exact answer, excuse by the whatever you can quickly do this and we get for N equal to the right answer. Similarly 3 of course there is no problem, 4 this is 4 plus 1 2.5 I should have at least 3 terms, so with 1 term I get 0, 2 terms I get 0.22 and this is 0.4, now 6, 6 plus 1, 7, 3.5 I would need at least 4 terms to get that, so 0.2857 is the exact answer, 1 doesn't give a good answer, 2 is this answer, 3 is this, 4 gives the exact answer, so you can quickly verify these things the information needed to check this table are here. Now what we, now we will return to the evaluation of the integrals, this is a mass matrix and using the first result I will be able to write in this form, upon transformation that is our first set of results that we got. Now how do we evaluate this integral is a question, now we propose to evaluate this integral using Gauss quadrature, so for that to do that correctly I should know the order of the function which is this integrand, that means how does N transpose N into determinant of J vary with respect to eta and Xi, if I know that then I will be able to select the correctly the number of integration points, otherwise I will again as you saw here if your integrand is X to the power of 6 and if you happen to take N equal to 3 you will not get the right answer, you have to take N equal to 4, because you know this rule you will be able to make this choice. Now we have to now work through the evaluation of this determinant of the Jacobian, so what we know right now Jacobian is given by this, Jacobian of the transformation X we are writing like this and Y we are writing like this, so given this I can now substitute for these derivatives in terms of NJs and using this expression I will be able to write J as product of two matrices, the first matrix is matrix contains derivatives of N1 and N2 with respect to Xi and eta and second matrix contains the coordinates of the nodes, okay, so all these quantities are known, so we can go ahead and perform the calculation, if you can it's quite straight forward to do this, if you do this I get J to be in this form where E1, E2, F1, F2 and F3 are explained here in terms of nodal coordinates, so moment I have this I can evaluate the determinant of J and it turns out to be this, and interestingly it turns out that the determinant of J in this case is a linear function of X and eta, it is precisely this that I wish to determine, then only I will know which order of integration I need to implement, we knew the order of N transpose N, because we know entries in N so we know what order that product would be, but we didn't know what was the order of this J is, I mean we didn't even know whether it is linear or not, first of all. So if we now look into this we see that it is a bi-coderatic or a bi-cubic function depending on which terms you are looking at, so using the rule that we are using a 2 by 2 array of integration points would be adequate to perform this integration, so mass matrix can be evaluated exactly by using Gauss quadrature with 2 by 2 Gauss quadrature point mesh, how about stiffness matrix, so I need here B transpose DB into the Jacobian, so again I need to perform some calculations here, so this D matrix has to operate now, so I need this information dou by dou psi is dou X by dou psi into dou by dou X plus dou Y by dou psi into dou by dou Y, etc., and this if you put in the matrix form I have this vector of gradient dou by dou psi and dou by dou eta, and this is the Jacobian matrix this, so I can either write this or the inverse of this I can express dou by dou X and dou by dou Y in terms of dou by dou psi and dou by dou eta through J inverse, now recall B was this for this problem, so we get after multiplication I get this and it turns out that the quantity that we are looking for is given in this form, okay, so elements of B transpose DB into determinant of J are ratios of bi-coderatic functions and linear functions, see there is a inverse here, so therefore if you look at individual terms they won't be polynomials, they will be ratios of polynomials, okay it turns out they are ratios of bi-coderatic functions and linear functions, so this would mean Gauss quadrature won't be able to evaluate elements of KE in an exact manner, so if you use Gauss quadrature to evaluate this we are introducing certain approximation, but Gauss quadrature is so convenient we will make some discussion on this subsequently, but right now the recommendation based on experience is used 2 by 2 Gauss quadrature, we will come to this slightly later once again, okay, so summary is mass matrix gets evaluated exactly and stiffness matrix there will be an approximation. Now this is, I have used now a four-noded quadratic element, now the question is if you are using discretization, if you are modeling any structural behavior how do you improve upon the accuracy? One way is stick to four-noded quadrilateral elements and use more of them, that means define the mesh, introduce more number of elements and do that, the other way is use a higher order polynomial to interpolate the field variables, if you want to do that it amounts to increasing the number of nodes in an element, so instead of having four-noded element we can have eight-noded element, okay, so then when I interpolate the interpolation functions will be of the higher order polynomials, they won't be what we saw for a four-noded element, because now there are additional nodes that come into our formulation. Suppose I take this, I will as before I will retain these three nodes and introduce additional nodes at the midpoint of each of these elements, again we will transform this to a square region and evaluate the elements of stiffness and mass matrices using Gauss quadrature, so this is a eight-noded quadrilateral element with 16 degrees of freedom. Now the field variable now needs to be represented in terms of the nodal values and there are eight nodal values for U, so if you go back here for a U1, U2, U3, U4, U5, U6, U7, U8, so I have to utilize all of that to interpolate U of X, Y, T you know to facilitate the development of this element. Now similarly V is represented like this, now using the logic that we already discussed we need not have to go into all the details every time, for we get these interpolation functions which have this coronacal delta property, for nodes 1 to 4 I get this is the interpolation function, for 5 and 7 these are the interpolation functions and 6 and 8 we get this. Now I would not like to get into all the details of the formulation of these integrants but if you indeed perform all that you will be able to notice that you need to use 3 by 3 Gauss integration points, with a 3 by 3 mesh mass matrix is evaluated exactly but the stiffness matrix continues to be evaluated approximately. Now how about this is a rectangular element, this discussion is not for this element, this discussion is for a rectangular element, for a quadrilateral element the same formulary works except now that the X and Y also need to be interpolated using the same shape functions, okay. So in this case it turns out we need to use 4 by 4 Gauss integration points in evaluating stiffness and mass matrices, and if you use 4 by 4 mesh the mass matrix is evaluated exactly, so the mesh size that you need to select should be such that the mass matrix elements of mass matrix matrix are evaluated exactly, see mass matrix elements turn out to be polynomials because N transpose N remains as polynomial, whereas here because of the Jacobian also of course would play a role, but here because of this you know various operations involved here the K matrix there will be difficulties as we saw while before. Now let us return to the shear wall and earth diameter example that we considered in the previous lecture, so the problem was there is a shear wall of thickness 0.2286 meters and these are the dimensions and the problem was to find, first find natural frequencies and normal modes. Now analytical estimates of natural frequencies using Timoshenko beam theory where obtained, pointed out these were their numbers, these are frequencies in hertz. Now we discuss this model in the previous lecture we use 64 elements with first order elements, that means we used four nodded quadrilateral elements, actually they were rectangular elements here and the degrees of freedom that the model had was 160. Now what I have done is we have retained the same mesh and same element configuration, but now the elements are modified to be second order elements, that means they are now eight nodded quadrilateral elements, so the degree of freedom increases to now 448 and the natural frequencies from 5.03 I come to 4.95 and these numbers change and you can see here they are reasonable, if this is believed to be correct because it's an analytical solution and the structure, structural geometry satisfies the requirements for the application of Timoshenko beam theory, so consequently we can test this result, so and this analytically derived, so we see that the first natural frequency is much better approximated than this, similarly the other frequencies also get represented better. So these are the mode shapes, first four mode shapes are shown here, so you can see that now the mode shapes are smoother because they are represented with greater detail within an element, this is the tenth mode, higher modes become representation of higher modes would not be so smooth as it is seen here, if you had used the first order element this mode would not be as smooth as it is seen here. Now this earth dam problem it was a triangular wedge model using plane strain approximation and these were the properties Poisson's ratio 0.45 this is the density and this is the angst modulus, with triangular elements with 42 elements and 42 degrees of freedom we got these natural frequencies and if you use the shear beam model and this is analytical estimate of the shear beam and natural frequency this is available in terms of Bessel's functions and then exact solution to this is available so based on that these numbers have been computed. Now this we have done in the last class, now what I have done now is this analysis is with first order triangular elements, now this analysis is second order triangular although I didn't discuss the formulation of second order quadrilateral element the logic is conceptually quite similar to what we briefly mentioned for quadrilateral elements, so if I use this same number of elements but the order has increased so degrees of freedom will increase, so the frequencies now change and the first frequency is if you have to believe this this is 1.25 this is 1.227 and this seems to be moving towards that, but mind you this is also an approximation to this behavior of this wedge, so the underlying assumptions of shear beam behavior must be satisfied by this structure, so it is difficult to pass a clear judgment on which of these two are exact. Now these are the mode shapes from a higher order element, this is the first mode and fifth mode, this matches with what we did with the Corsair element, I mean first order element, now in this illustration I am using a combination of triangle and quadrilateral elements, this is quadrilateral element, this here we are using triangular elements, so this is again first order elements, so there are 32 elements 56 degrees of freedom and we get natural frequencies like this, so this is just for illustration, and these are the mode shapes that we obtain from this model. Now what I do, I will repeat this exercise using the same mesh, instead of first order element I will use second order element, so this quadrilateral element here has 8 nodes and 8 degrees of freedom, the triangular elements here will have 6 nodes and 12 degrees of freedom, so number of elements remain the same, there are 210 degrees of freedom and I get different sets of estimates for natural frequency, so these are the mode shapes for the two first and the fifth mode. Now so far what we have discussed is basically free vibration problems, now if we have to consider first vibration response we need to worry about how the equivalent forces are computed, so let me just the basic rule for that is this, okay, so to explain this what we will do is we will consider a, suppose there is a edge 2, 3 in a triangular element and the loads will be in the plane of the element, the plane stress problem, suppose at any point here there is a load PX and load, sorry this load PY, now when we formulate the finite element model we also need to assemble the nodal forces, first of all we should evaluate the equivalent nodal forces, so what we do is we have this over the surface of the element delta U, delta V transpose PX, PY into DS, now we also have the approximation UV which is N into UV, so now we write delta WE for the nodal, in terms of the equivalent nodal forces these are the forces that are acting on the surface I want to represent them in terms of equivalent nodal forces, so we demand that equivalence of these two will give me the expression FE equal to, see I have to consider the edge 2, 3, I am assuming if you recall we are considering force on the edge 2, 3, so this has to be evaluated along that edge, so if I do this we will get FE for the example that we have considered as half length of 2, 3 into 0, 0, 0, PX, PY, PX, PY transpose, so what we are getting is at the two nodes half of the load comes here, half of the load goes there as per this formulation, so this is what the equivalent nodal forces, the logic for computing equivalent nodal forces, this can be done for quadrilateral elements and rectangular and quadrilateral elements as well, so at some point we will illustrate this with some examples, now in the subsequent lectures what we will do now, so far what we have done is we have dealt with a 3D beam element with two nodes and 6 DOF for node, now we have now considered 2D continuum, we have derived different elements, now the generalization for further analysis, now there are three things that we can do, this beam theory will get extended to plate and shell, this 2D continuum we can consider now three dimensional problems, that means when conditions of plane stress and plane strain are not satisfied, here another 2D model becomes possible that is axisymmetric solids, and then 3D solids, so what in the remaining next few lectures what we will do is we will develop the structural matrices for these problems and these problems, and at the end of it I will, at the end of these developments I will make a few remarks on choice of shape functions and their role on how the answers converge, now we will consider the question of how to compute equivalent nodal forces, we will quickly recall what we did for a beam element, so this is a beam element say 2 dimensional beam element deforming in its own plane, this is a distributed load F of X, T, so discretization of this element we had two nodes and two degrees of freedom per node U1, U2, U3, U4 were the four degrees of freedom for the element, the problem here is see the displacement field is represented in terms of nodal displacement and these interpolation functions, so how do we represent this distributed load in terms of equivalent nodal forces, so how to find P1, P2, P3, P4 in a way that these form an equivalent, there is an equivalence between this load and this set of loads, so what we considered was a virtual displacement, and that we represented in terms of virtual nodal displacements multiplied by the interpolation functions and the work done was equated, so the work done by virtual displacements on these set of forces and the virtual displacement on this, by this set of forces are equated, and by manipulating the terms we determined the nodal forces to be given by these expressions, so we are going to do something similar for even two-dimensional elements, so suppose PX and PY are the distributed loads along the edge 2, 3, now I want to find out what is the equivalent nodal forces because of these forces, so PX and PY are distributed along this edge 2, 3, so the virtual work is given by this expression, and the displacement field we are interpolating in terms of nodal degrees of freedom as shown here, now along this edge, so this virtual work can be written in terms of virtual displacement and applied forces as shown here, and this is our representation, so the virtual work which is this is also equal to nodal displacements into these forces, and after we manipulate this we get this expression, and we should note that this trial functions need to, here we need to evaluate this along the edge 2, 3, now this is the virtual nodal displacements into the nodal forces, and by equating these two we get the, by equating this and this we get the expression for equivalent nodal forces which is shown here, now here we have assumed that PX and PY are constant along the edge, otherwise there will be a quadrature that has to be implemented, the subscript L2, the quantity L2, 3 represent length of the edge 2, 3, this length, so if you know the nodal coordinates you can derive that length in a straightforward manner. Now we can make some remarks, first is the quadrature rule that we discussed was in the context of a scalar integral, this can be, the same formulation can also be extended to two dimensional integrals and one gets quadrature rule of this form. Now in this calculation of equivalent nodal forces we considered surface tractions, but if there are body forces within the element we can denote them by FX, FY transpose, suppose if that body force is constant within the element how we can find the equivalent nodal forces, but again we have to use the same formulation, I leave this as an exercise, this will be the equivalent nodal forces and here again we need to evaluate this integral and we need to examine the order of terms in the integrand and decide upon the quadrature rule, but this needs to be evaluated again using Gauss quadrature. Now a small exercise, so let's consider a four-noded linear rectangular element, now we consider the displacements in the discussion so far we considered this part, now to this I will now add two more terms alpha 1, 1-xi square, alpha 2, 1-eta square and similarly for V I will add this. Now you can see here that at the nodes which are xi equal to plus minus 1 and eta equal to plus minus 1 these terms are 0, so this is like a bubble function so they don't affect the behavior of U near the nodes, at the nodes and V is given by this similarly. So now again these trial functions are same as what we have used for a linear rectangular element, now these alpha 1, alpha 2 and alpha 3, alpha 4 are not nodal displacements now, they are generalized coordinates, so there will be now 12 degrees of freedom for the element, 4 nodal displacement U1, U2, U3, U4 and 2 generalized coordinates alpha 1 and alpha 2, this is 6 for U, similar 6 quantities for V. Now the exercise that is being suggested is use the energy expressions and obtain the 12 by 12 stiffness matrix and eliminate this alpha 1, 2, 3, 4 in terms of the nodal degrees of freedom using static condensation and obtain a 8 by 8 stiffness matrix, and also examine are the displacement fields continuous across the element boundaries, so this exercise can be carried out, so in using this element you may notice that we can continue to use the consistent mass matrix that was derived using only this approximation, kinetic energy can be characterized using only this, and whereas the displacement fields will be, the strain energy will be computed by an alternative representation. Now in the subsequent lectures now we will now consider problems of 3 dimensional elasticity, we have briefly discussed this earlier, so we will quickly recall, so stress, state of stress is now described in terms of 6 stress components, either it can be arranged as a matrix or a vector like this, these are the strain components, and for an isotropic linearly elastic material the stress and strain are related through this relation and this is the matrix D, and we have the expression for strain energy and kinetic energy, this is sigma transpose epsilon DV naught is the expression integral of that over volume element into half is a strain energy, and this is a kinetic energy. Now to compute strain energy we can now what we are doing is we are expressing sigma in terms of strain, and the strain subsequently we will express in terms of displacement as we have been doing, this is a strain displacement relation for small deformations, and substituting that we will get displacements as N into U E, where U E are the nodal displacements, and consequently epsilon is obtained in this form where B is equal to this matrix into N, so given that we will have the expression for strain energy and kinetic energy in terms of assumed displacement fields as given here, and we need to now carry out this integration over the volume elements. Now the volume elements in a 3 dimensional solid becomes more diverse, we can have a tetrahedron or a rectangular hexahedron or a pentahedron or an isoparametric hexahedron, and we can also have isoparametric versions of tetrahedron and pentahedron, and these can be first order elements or higher order elements, so there is a great diversity of elements, so we can examine some of them as we go along. Subsequently in 3 dimensional problems we can also consider solids of revolution, these are called axisymmetric problems, we will be doing that, and we will develop at least one element for this type of applications, and this would be followed by a discussion on a plate bending element, so this is a thin say lamina which carries now transverse loads in addition to the inline loads, so the action of this transverse loads on this lamina is to, the plate is to induce bending and about X and Y axis, and we need to formulate the correct expressions for strain energy and kinetic energy and develop the element, so we will take up these exercises in the following lectures, and may be some examples on force response analysis, so with this we will conclude this today's lecture.