 We have been discussing vibration analysis of planar frames, grids and 3D frames, so we will continue with that discussion, so in a planar frame the load resisting mechanism is essentially through bending and axial deformation, and the degrees of freedom in a typical beam element will have 2 axial deformations, 2 transverse displacements and 2 rotations. In a grid structure which again a planar structure the load acts normal to the plane of the structure, and the load resisting mechanism is through twisting and bending, so a typical grid element will have 2 twist degrees of freedom and 4 degrees of freedom corresponding to the bending action. In a 3-dimensional frame element, so this is a grid element where the field variables are theta and V, and we interpolate theta in terms of 2 linear interpolation functions and displacement V in terms of 4 cubic polynomials, so this is a grid element which is 2-noded element with 3 degrees of freedom per node, so this is the element stiffness and mass matrices for the grid element, and in a 3D element there will be axial deformation that is 1 and 7, twisting that is 4 and 10, and in plane bending and out of plane bending, so there will be therefore at every node there will be 6 degrees of freedom, so a typical beam element here will have a 2-noded beam element will have 12 degrees of freedom. Now in the last lecture we looked at the structure of the strain energy and kinetic energy for a 3D beam element, so the strain energy consists of energy due to axial deformation, due to twisting, due to bending, in plane bending, out of plane bending, where J is the torsional constant given in terms of you know psi which is a warping function. Kinetic energy there is kinetic energy due to axial deformation, due to twisting, due to the 2 bending actions, here IM bar is given by this integral. Now the field variables here are twisting axial deformation and 2 in plane bending out of plane bending, so twisting is represented in terms of linear interpolation functions, and similarly axial deformation is in terms of linear interpolation functions, V and W are in terms of cubic polynomials, so the Lagrangian in this case will be functions of U1, U2 to U12 is a 12 degree of freedom system, and application of this Lagrange's equation will lead to the element mass and stiffness matrix, and the element equilibrium equation as shown here. So the axial degrees of freedom are 1 and 7, and U of 0, T is U1, U of L, T is U7, and within the element we interpolate U of X, T using this function. Similarly theta is in terms of U4 and U10, and again linear interpolation function, V is in terms of U2, U6, U8, U12 using cubic polynomials, and W is in terms of U3, U5, U9, U11, and that's in terms of cubic polynomial. So if we now write the Lagrangian for this system and derive the element stiffness and mass matrix, it will consist of several components, so the component due to axial deformation will give us these 2 by 2 matrix, due to twisting it will give this 4 by 4 matrix, in plane bending this 4 by 4 matrix, and out of plane bending this 4 by 4 matrix. So all of that can be assembled, these are equivalent nodal forces in terms of the shape functions, here for example in the axial deformation case the equivalent nodal forces will be P1 and P7, and an axially distributed load, FX of X, T is represented in terms of equivalent nodal forces P1 and P7 using these relations. Similarly if there's a distributed torque about the longitudinal axis along the beam that is MX X, T the equivalent nodal torques will be P4 and P10, and that is given through this relation. So this is the usual nodal forces for bending action, this we have seen earlier. Now this is the configuration of the K matrix, A represents the axial deformations, T represents twisting, and BY and BZ are the in plane and out of plane bending respectively. So this is how the elements will be distributed in the stiffness matrix and also the mass matrix. So now if you populate these elements we get this 12 by 12 matrix, so you can see here AE by L minus AE by L, minus AE by L and AE by L is the contribution from axial deformation. Similarly GJ by L minus GJ by L, minus GJ by L, GJ by L is a contribution from twisting, and so in a similar fashion we can identify contributions from the two bending actions. Similarly the mass matrix can also be written following the same logic which is again 12 by 12 matrix, so this is displayed here. Now we have now obtained the 12 by 12 stiffness and mass matrix for a 3 dimensional beam element. Now the next question that we need to answer is what happens if the local coordinate axis for the beam element does not coincide with the global coordinate system which invariably happens in most applications. So to be able to do that we will now need to develop the rotation matrix in 3 dimension, so we will quickly run through this ideas here, so let us consider a global coordinate system capital XYZ and a local coordinate system lower case XYZ. Now this I cap, J cap and K cap are the unit vectors along XYZ, there is a global coordinate axis, this X cap, Y cap, Z cap are the unit vectors along the local coordinate system. A is a point in global coordinate system it has coordinates XYZ, upper case XYZ, in the local coordinate system it is lower case XYZ. Now therefore the position vector say OA bar I will call it as A bar can be written as XI cap, YJ cap plus ZK cap, similarly in the local coordinate system the position vector is given by this, if we want now component of the position vector along X axis we can consider the dot product between this position vector and this unit vector, so this is what we do, so that leads to this equation where cos XX is the angle between X cap and I cap, and similarly cos XY is a cosine of angle between X cap and J cap, right. So by using similar logic if we by taking dot product with Y cap and Z cap I get this relation, so therefore the relationship between the coordinates in local coordinate system and coordinates in global coordinate system is through this matrix, so we can deduce some properties of this T1 matrix, suppose if I denote the local coordinate vector XYZ as X bar and global coordinate vector as X bar and T1 as this matrix and the length of the vector OA as L, then this equation that we derived just now can be written in matrix form in this way. Now L square is a length of the vector which is given in local coordinates by this relation and in global coordinates through this relation, so this square of the length, so these two must be equal, okay, length is invariant with respect to the coordinate representation, so now for this X bar I will write use this equation and I get for X bar transpose I will write X bar transpose T1 transpose, so for this X bar I will write T1 X bar, so this must be equal to X bar transpose X bar by this relation, so if you observe now you see that T1 transpose T1 must be identity matrix, therefore T1 inverse is T1 transpose, or in other words T1 is an orthogonal matrix, now that is one important property of the so called rotation matrix, now if the two coordinate systems coincide then I will have this equation as X bar is equal to T1 X bar, that would mean lambda equal to 1 will be an eigenvalue of T1, okay, now this is another important property, now if lambda is an eigenvalue of T1 I have this eigenvalue problem T1 phi is equal to lambda phi, now if I pre-multiply by T1 transpose I will get T1 transpose T1 phi as lambda T1 transpose phi, now T1 transpose T1 is identity, so based on that I can reduce T1 transpose phi is 1 by lambda phi, so this would mean 1 by lambda is also an eigenvalue of T1, that means eigenvalues appear as reciprocals, now let's look at now the rules for coordinate transformation for the beam element, so the U1, U2, U3 denote the position of denote at this LF tent, so that is related to the position in global coordinate system let's call it as capital U1, U2, U3 through this matrix, similarly U4, U5, U6 is related to capital U4, capital U5 and U6 through this relation, on the right end 789 is related to the global 789 through this relation, similarly 10, 11, 12 is related to the global 10, 11, 12 through this relation, so this transformation matrix remains common for all this individual transformation, so if I now call vector U as U1, U2, U3 up to U12 and capital U as capital U1, U2, U3, U12 then this relation emerges from this set of 4 equations and I get this, where T bar is this matrix, this 0 is a 3 by 3 matrix, so this is a 12 cross 12 square matrix, now we can write U as TU therefore capital U will be T inverse U, but T inverse is T transpose, so it is T transpose U, now the expression for strain energy is half U transpose KU, now for where U this lowercase U is in the local coordinate system, in global coordinate system substituting for lowercase U in terms of this I get this equation, now by rearranging these terms you can show that K bar is T transpose KT, which is the stiffness matrix in the transformed coordinate system, similarly if we consider kinetic energy we can show that the mass matrix in global coordinate system is T transpose MT, now following the usual procedure for assembling and introduction of boundary condition the governing equation for the built up system can be obtained as shown here, this is typically this is how it will result, so you need to write the element equation of motion for every element and then take it to a global coordinate system, convert the node numbering to the global node numbering scheme and you follow the usual procedures and then impose boundary conditions partition the state vector into those components which are not known and those which are specified, and for the unspecified displacement vector this is a final equation of motion, so this we have done before a couple of times therefore we need not run into all the details once again. Now in implementing the computer in on a computer this transformation we need to evaluate elements of this matrix, this is cos XX cos XY etc., so we need to have an automated procedure for that, so that can be deduced as follows, so let us consider a 2-noded 3-dimensional beam element with 2 nodes named as I and J, so let me show the sketch this is a beam element this is I and J, this is a global coordinate system this is a local coordinate system, this XYZ and this is capital XYZ, now let XI, YI, ZI and XJ, YJ, ZJ be the coordinates of nodes I and J in the global coordinate system, so this I has coordinates XI, YI, ZI and similarly the other node has coordinate XJ, YJ, ZJ in the global coordinate system. Now we also select a point P in the local XY plane, right, and that point P we take it to have the coordinates XP, YP, ZP, okay, so you must understand that P is a point in the local XY plane but not lying on X axis, it should not lie on the X axis but it should lie on the local XY plane. Now our objective is now to get the direction cosines of this local coordinate axis, now as far as local X axis is concerned the direction cosines are given moment you know the coordinates of I and J you can reduce cos XX is this, cos XY is this and cos XZ, this is straightforward. Now how do you get the direction cosines of this local Z axis, so what we do now is we consider a plane formed by a vector along local X axis and another plane which contains the vector say the suppose this is O, so what we do is we can see this sketch, let me explain what are the entities on this graph, OI is a position vector of node I and this is given by this, OJ is position vector of the second node given by this, OP is a position vector of this point P which is given by this. Now to determine IZ we consider now two vectors namely a vector along the local X axis and this vector IP position, this IP, a vector along IP, these two vectors lie in the local XY plane, right, so any vector along the local Z axis must be normal to the plane containing these two vectors, so consequently I get IZ, vector IZ is IX cross IP, so this will help us to find coordinates direction cosines for IZ, moment I know the direction cosines of IX has been determined, IZ is now determined to find now IY what I do is I consider the cross product between IZ and IX and that should be equal to IY, so if we use these arguments we will be able to set up the required equations for the direction cosine, so let us quickly do that, so let IJKB unit vectors along global XY and Z axis, now we have IX cross IP is IZ, if I use this relation I get this, so now if you expand this determinant you can show that cos ZX is ZX by this capital Z, cos ZY is this, cos ZZ is this, where capital Z is this norm square root ZX square plus ZY square plus ZZ square, what are these ZX, ZY, ZZ they are given by the 2 by 2 determinants which are, which can be expanded like this, okay. So this helps us to find now cos ZX, cos ZY and cos ZZ, so cos XX, XY and XZ have been determined, so what remains now is the direction cosine of vector along Y axis, so let us consider Y be a vector along Y axis and the unit vectors X1 and Z1 along X and Z axis respectively, so clearly Y is a cross product of Z and X1, so now if I write that for Y I will write these components and the right hand side is the cross product, so all these cosines are known because we have already determined these quantities in the previous steps, so consequently I get cos YX, cos YY, cos YZ, in terms of components YX, YY, YZ and the norm capital Y modulus Y is square root YX square plus YY square plus YZ square, so this procedure can easily be implemented on a computer so what needs to be done is you have to select a point P in the local XY plane making sure that it doesn't lie on X axis and select P in the first quadrant of the local XY plane, that would help us in keeping track of science and getting the right answers. Now let me consider a couple of simple examples, so let's consider a 4-member 3-dimensional frame, so the coordinate system is XYZ shown here, so there is a member OD is along X axis, OA is along Y axis, OC is in the XZ plane, and OB is along Y axis, so the coordinates of these 4 points are shown here, we select the origin at O, and A is 3 meters, B is height is 4 meters, OD is 3.5 meters, C has coordinates 3.503. Now this is a properties I assume that the section is having a rectangular cross section although I have illustrated this for square, the formulation that we have is valid for general case. Now we name the degrees of freedom and we select, that we select has one element for every member, so there will be 6 degrees of freedom for each of these element at this common point, and they have to be shared therefore this structure in this particular model will have 6 degrees of freedom, 3 translations and 3 rotations at this point. Now we need to now form the stiffness matrix and mass matrix for this, so we have to be able to do the calculations for the global stiffness matrix, first we have to find the transformation matrix for each of the you know element, so I have given T1 for the position vector OA, OB, OC and OD using the procedure that we outlined you get this 4 matrices, then if by treating 1, 2, 3, 4, 5, 6 as degrees of freedom and assume imposing boundary conditions and so on and so forth, the reduced structural stiffness matrix which is 6 by 6 is obtained as shown here. A similar calculation is also made for the mass matrix and once we are ready with the reduced mass and stiffness matrix, matrices we can perform the eigenvalue analysis which I have done and I get this phi as the modal matrix as shown here and the natural frequencies we get as 61.76, 69.44, etc. these are in hertz, so this at least completes the free vibration problem and subsequently if there are any forces acting on the structure and if you have information about damping in the elements, etc. you could formulate the damping matrix, you can formulate the force vector and the complete equilibrium equation for force response analysis would subsequently be ready. So this is one some example of mode shape, the first mode shape I have shown, the red is the structure in its undeformed configuration, blue is the deformed configuration in first mode. Now this is a system with 6 degrees of freedom to understand how good this quality of this mode shape is, a more detailed model was also made and that detailed model had 366 degrees of freedom and the mode shapes, the first few mode shapes as per this model are shown here, so with respect to this, the discretization scheme that we have used this can be viewed as a, this is suddenly a more exact result, so you can see here that the first mode essentially consists of each beam behaving as if it's a proper cantilever. So the subsequent mode they are all in phase now, all the points are deflecting downwards and subsequent few modes each of the element will either be in phase or out of phase and those modes will appear and I have shown the first 4 modes. And mind you the coordinate system for these graphs are marked here. Now another example is again a structure with 4 members, so there is a rectangle AD CB and this is a O is a vortex and every A, B, C, D are joined to this through 4 members which has rectangular, each of which has rectangular cross section and we have taken the origin at D and these are the coordinates of the points O, A, B, C in meters that can be mentioned, and this is a coordinate system, so these are some of the parameters for the system. Now we again follow a simple discretization scheme for every beam, for every member I take one beam element, so the built up structure thus will have now again 6 degrees of freedom and all of them are located at this vortex, so 3 translations and 3 rotations. So now what we do, we form the element stiffness matrix and mass matrix for each of these elements and assemble, impose boundary conditions and get the reduced structural stiffness and mass matrices with which we will be able to perform the eigenvalue analysis, find natural frequencies and mode shapes which can subsequently be used for force response analysis. So here again for sake of completeness I have provided the transformation matrix for each of these members, position vectors for O, A, O, B, O, C, O, D with which you can construct the 12 by 12 transformation matrix for each of the members. Now after assembling I get this as a stiffness matrix and for this degrees of freedom being 1, 2, 3, 4, 5, 6. Similarly the mass matrix is a 6 by 6 matrix upon performing the free vibration analysis I get these as the natural frequencies, and since there is a symmetry in the structure the first two natural frequencies are pretty close, and this is the model vector. So this is a typical mode shape, the details are not very evident from this crude model but you can see that each one is again behaving like a propped cantilever, so to again to gain further insights into the nature of these mode shapes a more detailed model with 462 degrees of freedom was made, this is the first mode and as you can see all the four members are behaving like propped cantilevers, and the first few modes are essentially a combination of these four shapes being in phase or out of phase with respect to each of these members, so this is the mode shapes that we get. Now this T matrix for this problem for one of these members is shown here, now it turns out that we have shown that this individual 3 by 3 matrix is orthogonal, but it also turns out that this assembled transformation matrix is also orthogonal and I leave it as an exercise for you to figure out why this happens. Now the discussion that we have had till now we have treated beam as essentially as an Euler-Bernoulli beam, we have not paid attention to possible coupling between twisting and bending that happens in sections with unsymmetric, where unsymmetric bending is possible, and similarly coupling between axial deformation and bending that happens in curved beams, and there could be coupling between the two bending themselves, so there are many other considerations that we need to pay attention to if we want to complete the analysis of 3 dimensional frame. So what I will be doing is I will not run into each of this situation and explain what needs to be done, but on the other hand what we will do is we will quickly consider some of the complicating features, not all, and see what difference these considerations bring into our formulation, so I have selected a few issues, one is inclusion of rotary inertia that takes us to what is known as Rayleigh's beam theory, then presence of axial loads, then beams on elastic foundation, and beams which are rotating about an axis like turbine blades or ceiling fan blades which rotate about a vertical or horizontal axis, there are rotors in machines shafts which spin that I have not considered but that also is a special case that we need to think of, then deep beams and influence of shear deformation, see we have assumed the validity of Euler-Bernoulli beam theory where plane sections which are normal to the neutral axis before bending remain plane and normal to the neutral axis, so that virtually ignores, that means we are ignoring the shear deformation, contribution of shear to the transverse deformation. Now that assumption is reasonable if beam is not very deep, so if depth of the beam increases then these effects become important and how to include them in our formulation, then finally we ask the question is an exact dynamic analysis of 3D frames possible, we have asked this question in the context of a planar frames where we considered the dynamic stiffness matrix approach where the stiffness matrices were derived in terms of trigonometric functions, not piecewise polynomials, now the question I am asking is can a similar approach be developed for 3 dimensional frames, so let us see what are the answers for these some of these questions. So let's start with the so called Rayleigh's beam theory, now in the classical Euler-Bernoulli beam theory the displacement field that we assume U, so this is a coordinate system X is along the local centroidal axis of the beam, the assumption that plane sections that are normal to the neutral axis before bending remain plane and normal to the neutral axis after bending amounts to U being related to W through this relation, V is 0 that means we are not considering deformation along Y axis and the only deformation that we are assuming is W is W of X, T that means all cross sections we assume deformed by the same amount, now this is the assumed displacement field, from this now we can launch our calculations we can take first the strain displacement relations we can compute the strain, so the axial strain is dou U by dou X which is minus Z dou square W by dou X square, epsilon YY and epsilon ZZ are 0, how about the shearing strain? Dou U by dou Y plus dou V by dou X, so dou U by dou Y U is independent of Y that is 0, V is 0 therefore this is 0, gamma XZ is dou U by dou Z plus dou W by dou X, now dou U by dou Z is minus dou W by dou X dou Z and dou W by dou X is this, so these two terms cancel and there is no shearing strain, gamma YZ is again 0, so now the only stress component that we are having is sigma XX, okay which is this, all other stress components are 0 because all strain components except epsilon XX are 0 therefore corresponding stress components are also 0, now the strain energy therefore in terms of sigma XX and epsilon XX is given by this, and if we rearrange the terms I get integral over area Z square DA into other integral 0 to L E into dou square W dou X square DX, so this quantity inside the bracket is the moment of inertia area moment of inertia which is IZ, so we get the expression that we have been using all the while. Kinetic energy, this is where we have to pay attention now, in computing kinetic energy till in classical Euler Bernoulli beam theory we have not included the contribution from U dot square, so if you assume, if you look at the displacement field this is assumed displacement field, so when it comes to computing kinetic energy we are considering velocity only associated with W, we are not worried about U dot, right, now if you include that this will be the kinetic energy, so now U is Z dou W by dou X therefore U dot will be dou square W dou X dou T, so if you now substitute that into this expression I get this, so if you now rearrange the term the first term will be half integral 0 to L and this rho Z square DA I put inside this bracket and retain the other term outside, plus half again rho DA I keep inside the bracket and this is this, now this is the familiar mass per unit length which is M, now this is IZ, rho into IZ, so this is the contribution to the strain energy and we call this as rotary inertia, this is a new term in the mass, this will lead to new terms in the mass matrix which we have not included in the Euler-Bernoulli beam theory, so if beam is vibrating at high frequencies then this term becomes important, okay, so even for a beam which is not very deep as you go into the higher frequencies the cross sections will be deforming substantially and this will contribute significantly the rotary inertia, so this is a contribution made by the Rayleigh's beam theory. Now how about axial loads, so let's consider a beam which is carrying axial load, now because of these axial loads the deformed shape of the beam will be like this, so after deformation the beam will come somewhere here and this P into this delta will be the work done due to the axial load, so that we need to include in our calculation which we have not done because we are not simply considered that effect of P in the earlier calculation, so work done by the axial load is P delta, where delta is change in the length of the beam, how do you get that? So we can consider some section here, so this is DX, this is DW and this is DS, so from that triangle I can write DSS square root DX square plus DW square and consequently I can get this expression by taking out DX outside I get this and by using binomial theorem and retaining the first term I get this term, okay. Now so what is the length of the beam in the deformed configuration? I have to integrate this DS over the length, so that is given by this and if you see here I mean I am using the approximation I am ignoring the higher order terms so I get this, so this second term here is actually the change in the length of the beam, therefore delta is S bar minus L which is this, so the work done therefore is P by 2 into this, so now in formulating my expression for strain energy I need to include that, if I do that I get this expression, okay. Now kinetic energy if you include the effect of rotor inertia we will have this term and this is the kinetic energy due to the transverse vibrations, so this is what we can say is the expression for strain energy and kinetic energy for Rayleigh's beam which carries axial loads, okay, this is a new term, so now you can again you know W is the only field variable you can continue to use the Hermite polynomials and see derive your K and M matrix, I am not going to do that but it can be done, okay, that new matrix will have contributions in the stiffness matrix due to this axial load effect and contributions in the mass matrix due to contribution from rotor inertia, which are new, which is not there in our earlier formulation. We also have problems where beams rest on elastic foundation, so one of the typical model that we have for beams on elastic foundation is the so called Winkler's foundation model where the force offered along the length of the beam is proportional to the displacement, so now if I consider a beam resting on elastic foundation carrying axial loads, and I want to now develop a Rayleigh's beam theory for that, say. So this is the strain energy due to bending, this is the work done by the axial loads, this is now the new term which is the strain energy due to the elastic foundation. How about kinetic energy? Kinetic energy this is due to rotor inertia, this is due to transverse oscillation. Now the new term therefore is the contributions from the elastic foundation, this model can be improved upon, then there will be newer terms, okay, so this is one of the simple model for beams on elastic foundation. Again we can see that W is the only field variable, you can continue to use the cubic polynomials and derive the mass and stiffness matrix, so let me emphasize again in the stiffness matrix now there will be new terms due to the axial load effect and elastic foundations, and in the mass matrix there will be new terms due to rotor inertia, new with respect to the Euler-Bernoulli beam theory that we have discussed at length. Now let's consider a beam which is spinning about this axis as shown here with angular velocity of omega, now it is rotating the beam is spinning about the vertical axis here, about this axis. Now because of this there will be a centrifugal force which acts in the direction towards the center, I mean along this line and that can be computed by using this expression, okay, this gives centrifugal force at any point X along the beam I am ignoring the contribution from this eccentricity, the span is long compared to this eccentricity this can be ignored, so this is the axial load that will be present and that will be a function of this velocity of spinning, so now the strain energy will have energy due to bending and this is energy due to the spinning action, so this has to be now discretized, again W is, if you now consider kinetic energy this is the due to rotor inertia this is due to bending, so now our W is the only field variable we can go ahead and formulate the elements stiffness and mass matrix. Now we now consider transverse vibration of beams including shear deformation and rotor inertia, this is so called Timoshenko's beam theory. Now in this theory we assume that plane cross sections which are normal to the neutral axis before bending remain plane but do not remain orthogonal to the neutral axis, so there will be a shear deformation because of that, so this is the, that shear angle is shown as psi, so if you see the angle so the deformed thing will be something like that, so this angle between the neutral axis and this will not be 90 degrees there will be a shearing component. Now how do we formulate this problem? So we assume that actually we introduce new variable called psi that is at angle and U is given by minus Z into psi because the section is still plane it is not normal to the neutral axis but it is still plane therefore I can write this, V is 0 and W is W of X, T as before, so now we come, this is the assumed displacement form, now we consider the strains, epsilon XX is dou U by dou X I get this, epsilon YY and epsilon ZZ are 0, shearing strain gamma XY is 0, gamma XZ now I get this as the expression, similarly gamma YZ is 0, so this is a new term, so the stress components are therefore sigma XX is given by this and shearing strain sigma XZ is given by this, all other stress components are 0, therefore now in the expression for strain energy I will have contribution from sigma XX and sigma XZ and the associated strains epsilon XX and gamma XZ, so if you write that I will get this as the expression for the strain energy, the first one is due to the sigma XX, the second one is due to the shearing. Now we assume that the sigma XZ that is a shear stress according to the postulated displacement field is uniformly distributed across the depth, this of course is not true, we know shear stress is parabolically distributed across the depth for rectangular cross sections for example, so since we are assuming this to be constant we introduce a factor which multiplies the area of cross section and we call it as reduced section, okay, this is a constant, shear constant, so this is 5 by 6 for rectangular cross section and 1 by 1.175 for circular cross section, so using that assumption now the second term this term can now be simplified and I will get this expression for strain energy. Now there are two field variables W and psi, now the kinetic energy, kinetic energy is due to axial deformation and bending, so for axial deformation I use the postulated displacement field which is Z square psi dot square DA DX, and for the transverse vibration it is half this integral rho W dot square DA DX. Now we will rearrange these terms and I will put the terms rho Z square DA into a bracket integral over A and call that as IZ, and rho IZ is the, this is a contribution to the kinetic energy from rotary inertia, this is rho DA is mass per unit length which we have used earlier, this is contribution from inertia against translation. So now the Lagrangian will be in terms of psi and W, okay, now in principle we can develop the finite element model for this, okay, so this requires some careful consideration on what type of polynomials to be selected and other issues, so I am not going to get into the details, but this certain issue need to be examined before we select the right type of polynomials for representing W and psi, you must notice that the highest derivative that is present here is first order, say dou W by dou T and sorry dou W by dou X and dou psi by dou X are the terms having the highest derivative, so it looks like linear interpolation should work, but there will be some difficulties so you have to take care of that. So the final question that I asked was can we derive exact dynamic stiffness matrix for 3D frames, that means exact frequency response functions where there is no error due to model truncation and there is no error due to the use of piecewise polynomial interpolation functions, the word exact is within the framework of a given beam theory, okay, if you are using Euler-Bernoulli beam theory it is assumed I mean you have to accept that that theory itself is approximate, but within the framework of that theory we are getting exact solutions. Now the question we are asking is how to get this dynamic stiffness matrix for 3D frames, the dynamic stiffness matrix without using polynomial interpolation can also be derived for a 3 dimensional beam element. The resulting formulations again lead to exact solution for system FRFs for grid and 3D frame structures, this we have shown for plane frames, the simple observation that we are making now is that it can be extended to grids and 3D frames. So we can quickly run through the arguments, for axially deforming bar we have derived the dynamic stiffness matrix which was this, where lambda 1 was a parameter defined in terms of mass, axial rigidity and structural and viscous damping terms. Analogously we can also derive the dynamic stiffness matrix for a bar undergoing torsion, harmonic torques responding harmonically in steady state, so this is the dynamic stiffness matrix. The lambda 2 now will be different from this, instead of M I have IM bar instead of AE I have GJ and instead of H2 C2 H1 C1 I will have H2 theta, C2 theta H1 theta and C1 theta because these are the damping terms against twisting, whereas these are damping terms against axial deformation, they need not be the same. Bending we have done this, now this is in plane bending the degrees of freedom are 2, 6, 8, 12, so we have derived this dynamic stiffness matrix and this has a parameter lambda, I call it a lambda 3 here, this again is now function of mass per unit length and flexural rigidity EIY, similarly the out of plane bending we can formulate the dynamic stiffness matrix as we have done for a single beam, that brings in a new element, a new parameter lambda 4, which is again given in terms of mass per unit length, flexural rigidity and damping parameter. So all this can be assembled into this format which we have done for static stiffness matrix, so the structure of the dynamic stiffness matrix would essentially be the same as that we got for static stiffness matrix, so that would mean there will be 4 terms corresponding to axial deformation, 4 terms corresponding to twisting and 16 terms for bending in Y, about Y and bending, 16 terms for bending about Z, so this will be the structure of the stiffness matrix. So for a built up structure such as the 2 examples that we considered, for each of the element we again formulate the dynamic stiffness matrix, the coordinate transformation rules are exactly the same as has been used for static problems, and the reduced equation now will be an algebraic equation, so that can be inverted for every omega to get the required frequency response functions. So in principle we will be able to obtain the frequency response functions without resorting to modal truncation and without having to use piecewise polynomials for interpolating the field variables within an element, so this type of solutions can serve as benchmark solutions to validate approximate methods, and these are quite flexible in terms of their ability to handle different damping models, so if you are doing a damping treatment for individual members etc., this formulation will give you a better hold on the problem. So we will not pursue this, this is somewhat a special topic but it is useful to note that such solutions exist. There are other considerations which I have not talked about, there will be rotating shafts as in rotating machinery they will undergo vibrations, and then there will be problem of bending torsion coupling, for example in a channel cross beam with a channel cross section bending in one direction is typically accompanied by twisting, and similarly in curved beams bending and axial deformations will be coupled, and in certain cross section there will be coupling between bending about the two orthogonal axis, and there would be beams with pre-twist as in certain machine elements, turbine blades for example there will be a predetermined pre-twist, and how do you include that in your formulation? Then beams could be laminated, you can use higher order theories, there could be piezo materials inside embedded in the beam so on and so forth, so this is a list that can continue for various, I mean continue to address various issues that we have not addressed, but it is not crucial for us to go into every bit of these details, but my idea of pointing this out at this stage is these are issues that we are not covering, so if you are interested the discussion that we have made till now hopefully you would serve as a starting point for you to pursue these topics. Now we have reached now a stage in the course where we are able to formulate the equation of motion for skeletal structures, three-dimensional skeletal structure, now what we will do is now we will switch back to some issues about numerical solutions, how to solve the equations of motion, the approximations that we have introduced till now has been only in space, time has remained as a continuous variable, so now to solve the governing equations we have to integrate this set of coupled ordinary differential equation, so at the end of modeling that we have been describing the final equation of motion will be of this form, MU double dot plus CU for sake of generality I have introduced certain nonlinear terms which has not emerged from analysis that we have done till now, but you could imagine that if you include nonlinear strain displacement relations or nonlinear stress strain relations then the equation would no longer be linear and it may have terms like this, F of T is a generalized nodal forces forcing vectors and these are the initial condition. Now this equation constitutes a set of semi discretized system of coupled second order ordinary differential equation, semi discretized means in space we have accomplished discretization but in time it is still, time is still continuous it has not yet been discretized that is why we call it a semi discretized, that is these equations have been obtained after discretizing the spatial variable, this set of equation constitutes a set of initial value problem, that means that the upon integrating these equations there will be arbitrary constants and in order to determine those arbitrary constants we know the state of the system at a same value of the independent variable namely T equal to 0, therefore they are called initial value problem. So now the next task on our hand is to develop numerical methods to integrate these equations of motion, so what we will do is in the next class we will consider this question and try to develop numerical schemes which would enable us to obtain time histories of displacement and velocities, and subsequently consider questions on how to evaluate stresses and other quantities of interest, so with this we will conclude the lecture at this point.