 We have been talking about the development of finite element method for analyzing planar structures, so we will continue with the discussion. In the last class we outlined the basic framework for finite element modeling, so if omega is a domain that is of interest to us, we partition the domain into a set of subdomains called elements, and the union of these elements approximately equals the domain of interest, and with this each domain is called an element, and what we do is we approximate the field variable that is U of X1, X2, X3, T in terms of values of the field variables at the nodes, at a set of N nodes they are these dots shown here, and within an element we approximate the field variable in terms of the values at the nodes and by using interpolation, this NIE of X tilde are the interpolation functions. So the collection of these elements which approximates the domain of interest is known as the mesh, so this approach basically leads to an approximate numerical method to obtain solutions of problems posed as either in terms of a partial differential equation or in terms of a variational principle. Now in the last class we considered the element behavior where the element was modeled as an axially vibrating rod element, so it's a line element, these are the parameters of the element, axial rigidity, mass per unit length, length and damping, and we assume that the cross section is symmetric specifically we took that it is rectangular, and we identified two nodes X equal to 0 and X equal to L, and the value of the field variable at X equal to X1 was U1 of T, and the value of the field variable at X equal to L was U2 of T, and within this rod that the displacement field at any point X is expressed in terms of the nodal values U1 of T and U2 of T, and these two functions 1 minus X by L and X by L are the interpolation functions. We showed that using Lagrange's equation that means assuming this is a displacement field we can substitute this into the expression for kinetic energy and the potential energy and write the Lagrangian and apply the Lagrange's equation we get the element level equation of motion as shown here, and this U1 double dot and U2 double dot are the acceleration vector is the acceleration vector and U1, U2 this is a displacement vector, this M which is this is the element mass matrix, this is the element stiffness matrix, I also pointed out that we have used the same interpolation functions in computing elements of mass matrix as well as elements of stiffness matrix therefore this mass matrix is called the consistent mass matrix. Now we will continue with this discussion and we will now consider the behavior of an element which is modeled as an Euler-Bernoulli beam, so here we have again this is a line element this is X and Y are the coordinate axes and the displacement field along Y is denoted by V, V is a function of X and T, and we make some of the standard assumptions that are made in Euler-Bernoulli beam theory, we are considering bending in the XY plane and we assume that X axis coincides with the central axis of the beam, and we also take in this study that there is no coupling between bending and torsion, a possibility of including this can be considered later, but right now we will assume that there is no unsymmetric bending in the problem. So based on the standard assumption of Euler-Bernoulli beam hypothesis that for example plane sections which are normal to the neutral axis before bending remain plane and normal to the neutral axis after bending, we can postulate the form of the displacement field, so this is U in terms of V, and using this we can compute the strains and that can subsequently be used in the expression for strain energy, and we get after a few simplifications the well-known expression for the strain energy due to bending. So IZ is the moment of inertia, area moment of inertia about the Z axis this is given here. Kinetic energy is mass into velocity square integral over the domain half of that and if we use now rho A, integral rho A over A that is mass per unit length we get the kinetic energy in terms of mass per unit length and the velocity field. Now based on this we can now construct the Lagrangian which is given in this expression, so this is the kinetic energy and this is the strain energy. Now to obtain a finite element discretization we will now take that the element is made up of two nodes and at each node the translation and the gradient that is a slope are taken to be the nodal degrees of freedom, this is U1 and U2 of T are the values of, U1 of T is the field variable displacement at X equal to 0, this is dou V by dou X at X equal to 0, similarly U3 of T is V of X, T at X equal to L and this U4 of T is dou V by dou X at X equal to L. Now this element is a two-noded element with two degrees of freedom per node and it has four degrees of freedom. Now the displacement field within the element as before we are going to now express in terms of the nodal degrees of freedom U1 of T, U2 of T, U3 of T and U4 of T as shown here. Here phi 1, phi 2, phi 3, phi 4 are the interpolation functions which we need to now select. Now this UI of T for I equal to 1, 2, 3, 4 are the generalized coordinates and this phi of X are the trial functions. Now let us impose the boundary conditions that we already know at X equal to 0 and X equal to L, at X equal to 0 we know that V of 0, T is U1 of T. Now that would mean V of 0, T is now if I substitute in this expression it is U1 of T, phi 1 of 0, plus U2 of T, phi 2 of 0, etc. Now I know that this must be equal to U1 of T, so that can be achieved by inspecting if I now take phi 1 of 0 to be 1 and phi 2 of 0, phi 3 of 0 and phi 4 of 0 to be 0 I will be able to satisfy this requirement. Similarly if we now consider dou V by dou X at X equal to 0 I know that it is U2 of T. Now let's again use that V prime of 0, T is in terms of this, in terms of the assumed representation it is given by this equation. So again in order to satisfy this we can select phi 1 prime of 0 to be 0, phi 2 prime of 0 to be 1 and this phi 3 prime of 0 and phi 4 prime of 0 to be 0, so this will ensure that this condition is met. Similarly I can impose condition at X equal to L where V of L, T is U3 of T and here we will require phi 1 of L to be 0, phi 2 of L to be 0, phi 3 of L to be 1 and phi 4 of 0 to be 0. Finally dou V by dou X at X equal to L this is U4 of T and to achieve this we have to take this conditions on phi's that is phi 1 prime of L is 0, this is 0, this is 0 but this has to be 1. So this choice if we now select phi of X so that this set of 4 conditions which represent 16 conditions in total if we satisfy that we would ensure that we have satisfied the requirements on displacements at X equal to 0 and X equal to L. Now we can to now make the choice for phi of X we can now consider the equilibrium of the beam under the action of support displacements UI of T I equal to 1, 2, 3, 4. So for phi 1 for example we have 4 boundary conditions, so now if I write this equation for phi 1 we get phi 1 of X to be a cubic polynomial and for phi 1 I have 4 conditions that we have derived, now if I impose that I will get phi 1 of X to be this function, using the same logic that means as you make a cubic polynomial with 4 constants A, B, C, D we can derive phi 2 of X, phi 3 of X and phi 4 of X, okay. Now how do they look like? So phi 1 of X is this blue line so at X equal to 0 it is 1, its slope is 0 and at X equal to L the value is 0 its slope is also 0, similarly I have phi 2 of X the slope here is 1 and the value is 0 and here again the value and the slope are 0, so we define phi 3 and phi 4 in the same way. So these are the interpolation functions that we use to represent the beam behavior. Now I have now the approximation for the displacement field in terms of nodal values U1, U2, U3, U4 and these 4 interpolation functions, so now I can substitute this into the expression for kinetic energy and the potential energy, suppose now if I substitute into the kinetic energy I have this is the expression for kinetic energy, so this representation leads to this and noting that these phi I of X are all polynomials you will easily see that in this case we will be able to integrate, carry out these integrations in closed form, so the expression for kinetic energy in terms of the nodal degrees of freedom we obtain as shown here, so you have to carry out these integrations which are straight forward. Similarly the expression for strain energy is given by this and here I substitute for dou square V by dou X square which is UI of T phi I double prime of X summed over I equal to 1 to 4, so this again phi I double prime of X will again be polynomials and therefore carrying out this integration is a straight forward exercise and I will get this as an expression for the strain energy. So now my Lagrangian is ready T-V in terms of nodal degrees of freedom, there are 4 nodal degrees of freedom, so this is a system with 4 degrees of freedom and 4 generalized coordinates, so I need to apply the Lagrange equation on these 4 variables, so if I do that I get the equation in the form MU double dot plus KU equal to P, where U is a 4 cross 1 vector, U is 4 cross 1, U double dot is 4 cross 1, so this is the vector of nodal displacement, this is a vector of nodal acceleration, M which is a 4 by 4 symmetric matrix is the element consistent mass matrix, K which is this 4 by 4 symmetric element stiffness matrix, so the M is a element consistent mass matrix, it is symmetric, K is element stiffness matrix which is symmetric, so the element that we are talking about therefore is a 2-noded beam element with 2 degrees of freedom per node, and K and M we see that they are non-diagonal, and we can express the expression for strain energy and kinetic energy can be obtained in terms of nodal displacements and nodal velocities as shown here, this can be verified. So we are now being able to derive the mass and stiffness matrix for an Euler-Bernoulli beam element. Now again let us go back to the nomenclature, I call this mass matrix as consistent mass matrix, this is there is another way of forming mass matrix which is known as lumped mass matrix, for example if I, if this is a beam element and if half of this mass is lumped here and this half is lumped here I get a mass matrix which is as shown here, the half the mass is lumped here, this is one and half the mass is lumped here, this is one, so this obviously doesn't take into account various cross coupling terms that are present in the consistent mass matrix. Now we will consider some simple examples so that we understand how to use these matrices, suppose if I have a property cantilever, so I model this element as this beam using only one element, therefore there is a 2-noded beam element with 2 degrees of freedom per node, now the imposition of the boundary conditions for this problem results in this statement U1 of T is 0 because this displacement is 0, U2 of T is 0 because the slope is 0 it is a clamped end and similarly U3 of T which is the displacement here which is 0, P4 of T which is the bending moment here which is 0, so now if I write the element level equation of motion I get now for the acceleration vector the first 3 terms would be 0 and similarly for the displacement vector the first 3 terms will be 0, this U4 is the unknown associated with displacement which are 0 at the boundary there will be a stress reaction, for example there will be shear force and bending moment here and there will be shear force here and there this R1, R2, R3, so this equation therefore is a set of 4 equations, 3 for the reactions and 1 for the unknown displacement, so now we solve for the equation for the unknown displacement I get this, so you have to use the equation corresponding to the last row of this equation, so you get this single equation. Now this is a single degree freedom system, in this approximation the beam is approximated as a single degree of freedom system and the generalized coordinate is U4 of T, once I know U4 of T I will be able to write the expression for the beam displacement, how do we do that? That is all, 5, 4 of X is the polynomial, cubic polynomial that we have already derived, so this is a displacement field for all values of X and this U4 of T itself is governed by this differential equation, so from this you can deduce the natural frequency of the system, this is a single degree freedom system, square root K by M is the natural frequency, if we do that we get this as the answer 20.49 square root EI divided by ML to the power of 4. Now an exact solution to the natural frequency of this beam is available and if we, this is reproduced here and by comparing these two we see that in this approximation the first natural frequency computation has an error of about 33%. Now a similar exercise we can do for a simply supported beam, again this beam I approximate using only one element, so here you see that U1 is 0 because this is a simply supported end, U3 is 0 and P2 of T, what is P2 of T? The bending moment here which is 0 and P4 of T which is a bending moment that is 0, so what are the unknowns in this problem? U2 and U4 and the reactions shear force at X equal to 0 and X equal to L which is P3 and P1. So now if you write the element level equation of motion we get this is the equilibrium equation for the element, now for the beam now because we have used one element approximation so this is U1 and U3 are 0, U2 and U4 are nonzero, so whereas R1 and R3 are the unknown reactions which we have to compute, so again this is a set of four equations, two for the unknown displacements U2 and U4 and two for the unknown reactions R1 and R3, now if you extract from this the equation for U2 and U4 I get a 2 degree of freedom approximation to the simply supported beam structure. The remaining two equations govern, you know they are the equations for the unknown reaction, so though if you want to find the reactions you have to first solve the equation for the displacements and substitute here. Now let's now look at this 2 degree of freedom system and if we are interested in finding natural frequencies we have to solve the eigenvalue problem, this is the statement of the eigenvalue problem that means we assume that the two generalized coordinates oscillate harmonically at the same frequency omega and for that to happen this condition must be satisfied and if we now introduce a parameter lambda that means I will divide by EI by L cube and call this parameter as lambda the eigenvalues now expressed in terms of lambda turn out to be 2 by 7 and 6 from which I will compute the first and the second natural frequency, it's a 2 degree freedom system therefore I will get an approximation to first two natural frequencies, so these turn out to be 10.95 into this factor and 50.20 into this factor. Now again this problem is well known to have the exact solution and the exact solution is for the first natural frequency instead of 10.95 I have a factor 9.87 instead of 50.20 I have a factor 39.48 and that would mean in this model the errors in the first two natural, calculation of first two natural frequencies is about 11% and 27%. How do we improve upon this approximations? We have to introduce more nodes and increase the size of the problem, so we will see such issues later. Now so far what we have done is we have studied an actually vibrating bar and a flexurally vibrating beam, a single span, a single element behavior has been studied. Now if you consider a planar structure like this carrying P of T, if you carefully you know see each member of this structure not only deforms axially but also it flexes simultaneously, that means the bending and axial vibrations coexist in each of these elements, or in other words if you compute strain energy or kinetic energy in the each of these elements there will be contribution from axial deformation as well as contribution from flexion, so how do we treat this problem? So we will now address that issue, this leads us to the discussion of what is known as 2D beam element, so 2D beam element with you know shown here again we assume that there are two nodes, now the nodal degrees of freedom for axially vibrating rod it was only the axial deformation, now for only bending it was the translation U2 and the rotation U3, now since both axial deformation and bending coexist at every node I should now consider U1, U2 and U3, one axial translation and one transverse displacement and one rotation, so I have U1, U2, U3 at one node and similarly U4, U5, U6 at the other node, so therefore this beam element has two, it is a two-noded element with 3 degrees of freedom per node, it is a 6 degree of freedom element. What are the parameters of this model? The flexural rigidity EI, axial rigidity AE, M is mass per unit length, L is length and C is a damping parameter, we again assume that flexure is according to Euler Bernoulli beam theory, so the field variables in this case are U of X, T which is axial deformation and V of X, T which is the transverse displacement, so these two are the field variables, these are the nodal values of these field variables, for example U1 and U2, U4 are the values of U of X, T at X equal to 0 and at X equal to L, similarly U2, U3, U5, U6 are nodal values associated with V of X, T and dou V by dou X at X equal to 0 and at X equal to L. Now we need to now again interpolate the value of these field variables V and U in terms of the nodal values, how do we do that? We will see that, now the expression for kinetic energy and potential energy, if you now consider the strain energy has now contribution from axial deformation and contribution from flexure, so we have to add them, similarly the kinetic energy has contribution from axial deformation and transverse deformation, so these are the expressions for the strain energy and the kinetic energy. Now the axial displacement field is again represented like this, U1 and U4 are the nodal degrees of freedom and phi 1 and this should be phi 4 are the interpolation functions and phi 2, phi 3, phi 5, phi 6 are the interpolation functions associated with bending, so phi 1 is 1-X by L and phi 4 of X is X by L, these are the linear interpolation functions associated with axial deformation and these are the cubic polynomials associated with bending. So now I can form the Lagrangian which is V-T, now this will be function of 6 generalized coordinates, U1, U2 up to U6, so this is a 6 degree of freedom model, so I can write now the Euler Lagrange equation and for the 6 variables I equal to 1 to 6 and if I do that again I get a matrix form of equation of motion for the element, and this M for this case turns out to be this 6 by 6 matrix, and similarly the stiffness matrix is in terms of this 6 by 6 matrix, it is again symmetric and you can see that the contribution from axial deformation is contained here, and if you partition and remove these rows and columns what remains is the consistent mass matrix for flexure, similarly here if you look at this 2 by 2 matrix it is the stiffness matrix corresponding to the axial deformation, this R is the radius of gyration that is square root I by A that is how we have written this, and the remaining part suppose if you delete these rows and these columns we get the 4 by 4 matrix which is a stiffness matrix of the Euler-Bornley beam that we have derived. Now so if we now return to this example one complicating issue namely that in each of these elements we have coexistence of bending and axial deformation we are able to take care of, so now if you want to discretize this beam and study this you can use for each of the element that you choose a representation of this kind. Now there is another problem here there are few more problems that we need to address one by one, the next problem I would like to consider is now if you look at these members associated with each of these members I can think of a local coordinate system, and for the entire structure I can have a global coordinate system, for example this could be the global coordinate system whereas for this element the local coordinate system will be this. Now how do we now take care of the fact that the local coordinate system for this element does not agree with the global coordinate system, at the end of the day we should be able to describe the behavior of this structure in a common coordinate system, so now let us address that issue, so that takes us to the question of transformation from local to global coordinate system, so this is you know we have to invoke the standard results from transformation of coordinates as applied to vectors, so suppose you have a coordinate system X and Y called this as global coordinate system, and a local coordinate system lower case X and Y and this orientation is a gamma, then a point here viewed from this X Y coordinate system will have coordinates X and Y, if you view from the local this local coordinate system it has coordinates X and Y, now we can ask the question how this capital X and Y and lower case X and Y are related, so a simple consideration of the you know geometry here leads us to the relation lower X is X cos gamma plus Y sine gamma, and Y is minus X sine gamma plus Y cos gamma, so in the matrix form I can write this as X Y in local coordinate system is related to X Y in global coordinate system through this transformation matrix, and now if we use notation M and N for cos gamma and sine gamma I get this coordinate transformation matrix, suppose if I call it as C, the C is M, N, minus N and M, as is well known this C matrix is orthogonal that means C transpose is C inverse, now we will basically put to use this result to achieve what we want, so now let us consider a beam element which is like this initially we will view the element in the local coordinate system, so U1, U2, U3, U4, U5, U6 are the degrees of freedom with respect to the local coordinate system, in the global coordinate system they are U1 bar, U2 bar, U3 bar, U4 bar, U5 bar, U6 bar, so now if you take the location of this node the coordinates in local coordinate system is U1, U2, in global coordinate system it is U1 bar, U2 bar, so now you go back and use this result for this I get now I am going from local to the global so this matrix is replaced by transpose, so this is the position coordinate of this node, this is how the transformation works, similarly if you come to this node the position coordinate is U4, U5, and it goes to U4, U5 bar through this relation, now the rotation about Z axis remains the same because Z axis has not undergone any transformation, so U3 bar remains as U3, U6 bar remains as U6 therefore this transformation matrix is an identity matrix. Now we can collate all these results and put in a single form the 6 degrees of freedom namely U1, U2, U3, U4, U5, U6 in local coordinate system is related to the corresponding values in global coordinate system through this transformation matrix, so this is 6 cross 1, this is also 6 cross 1, therefore this transformation matrix is 6 cross 6, so we can write this in a compact form, I write this as U bar as T0 transpose U, where T0 transpose itself is written in this form we notice that these matrix this is L, this is L that repeats and this 0 is the 3 by 3 0 matrix this, so this is L, so I can write for U as T0 U bar, this is local this is global. Now what we need to do is use this transformation in our expression for the kinetic energy and potential energy, so in local coordinate system the kinetic energy is given by this, the coordinate transformation doesn't change the energy, the energy is a scalar, but the representation for mass matrix changes, so that is what we are looking at, so I know that U is T0, U therefore I can write U dot transpose as U bar dot transpose T0 transpose MT0 this, so this quantity which is T0 transpose MT I call it as M bar, so this is a mass matrix in the global coordinate system for the element that is given by T0 transpose MT, and the kinetic energy now I have in terms of the velocities expressed in the global coordinate system, similar exercise we can now repeat on expression for strain energy, so I have this as half U transpose KU, for U transpose now I will write, for U I will write at T0 U bar and U transpose will be this, and again I introduce a K bar which is T0 transpose K T0, so this is the element stiffness matrix in the global coordinate system, and this is an expression for strain energy in the global coordinate system, so I have now therefore K bar and M bar which are the structural matrices, element stiffness matrix in the transform coordinate system, M bar is the element mass matrix in the transform coordinate system, so we have taken care of this complication that the global and local coordinate systems could differ for different elements. Now we will now move to the next complication that we need to deal with, now I will list what are the complication we need to deal with at this stage, we need to now take care of the issue that this structure is made up of different elements which are oriented in a different direction, so if I now partition this into elements I will probably take this as one element, this as at least as one element, this as one element, this as one element, we could have more elements within this, that's a different issue, so the question is how to analyze a structure which is built up of different elements, that is the first question, the next question we have to deal with is how to deal with external forces which I have not yet tackled, and how to impose the boundary conditions systematically, and if there are any support motions as in case of earthquake engineering problems how do they enter our formulation, and finally I have not talked about damping, how do we handle damping in this system, so we have to answer all these questions so I am taking up one by one, so now the next question I wish to answer is how to analyze built up structures, or in other words if I now assume that this is one element this is another element this is one element and so on and so forth, for each of these elements I can derive a 6 by 6 mass and stiffness matrix, but how do I now get the mass and stiffness matrix for the built up structure by assembling these element level structural matrices, so how do we do that, to illustrate that we can consider a simple problem suppose there is one beam which has been divided into two elements, say S equal to 1 and S equal to 2, now let this be a Euler Bernoulli beam there is no axial deformation in this illustration, so how many degrees of freedom this structure has it has 6 degrees of freedom, this is 1, 2, 3, 4 and 5, 6, now each element has 4 degrees of freedom, so now if I take this element separately I will number these degrees of freedom as 1, 2, 3, 4, if I take the second element separately this numbering as per our scheme will be again 1, 2, 3, 4, but we should somehow impose the now the condition that what is 1 here must be 3 here, what is 2 here must be 4 here that would ensure compatibility of displacements, how do we do that, so the basic idea is the displacement at shared nodes must match, so let P be the number of elements which was 2 in the illustration that I showed just now, let N be the total degrees of freedom which was 6 in the illustration, now I will denote the degrees of freedom as UJ bar J running from 1 to N, now the nodal displacement vector is a N cross 1 in global coordinate system for this structure, now this is UJ bar J equal to 1 to N, but for each of the elements as I have shown here I have this lower case U bar of S which are the nodal displacement vector in global coordinate system for the S element, S equal to 1 for example I had 1, 2, 3, 4 and S equal to 2 again I have 1, 2, 3, 4. Now we would like to write the nodal degrees of freedom at the element level in terms of the nodal degrees of freedom at the structure level, so what we do is we introduce a transformation U bar of S is written as AS into U bar, now U bar of S is some NS cross 1, in this case NS is 4, each beam element has 4 degrees of freedom it is 4 cross 1, whereas this U bar is 6 cross 1, therefore this A matrix should be NS cross N, so number of rows in A matrix will be equal to the number of degrees of freedom at the element level, and number of columns in A matrix will be equal to the total degrees of freedom at the structure level, so now for the simple example that we used for this example if you now consider U1, U2, U3, U4 for the first element what is 1, 2, 3, 4 is same as 1, 2, 3, 4 for the global structure by that I mean this 1, 2, 3, 4 and this 1, 2, 3, 4 and this 1, 2, 3, 4 are the same there is no difference, so the A matrix in this case will have an identity matrix here and then zeros, for the second element if you go back what is 1 is actually 3 for the global structure and what is 2 is 4, and similarly what are 3 and 4 are respectively 5 and 6, so the A matrix now will become as shown here, so U1 bar is the third element U3 bar, U2 bar is the fourth element which is a U4 bar and so on and so on. So with this representation now we have the expression for kinetic energy and strain energy in terms of the local node numbering scheme now we will go to the built up structure, now the strain, the kinetic energy for example in the built up structure consists of contribution to the total kinetic energy from each of the elements, so if there are P elements this is some of S equal to 1 to P TS of T, now TS of T itself is given in terms of the local degrees of freedom, the element level degrees of freedom as shown here this is what we have derived now, these are in the global coordinate system. Now I will introduce the transformation U bar is AS capital U, so if I substitute this I get this expression, now this summation is over the elements and you can see that the quantities inside the braces here and here are independent of S, so they can be pulled out and the summation will be only on these terms, this summed up quantity which is written as this is the structure mass matrix in the global coordinate system, okay. So this has, this receives contributions from each of the element, that is the element level mass matrices combined together according to this rule to get, to lead to the global level mass matrix. Similarly the strain energy is again the total strain energy is sum of contribution from strain energy within each elements and this is again following the same logic it is given in terms of the element level degrees of freedom and in terms of global degrees of freedom through this transformation like this, and again we can pull out these quantities in the brace outside the summation and I get this and this summation is actually the assembling procedure leading to the global structure stiffness matrix. So this is summed over S equal to 1 to P this, so we have now therefore the expressions ready for global mass matrix and global stiffness matrix. Now we can write notionally the equation of motion that is a global equation of motion, the right hand side we will call it as a force vector which we have not yet you know describe what it should be, hmm. So now let us come to that description but before that we can you know illustrate what exactly we have done by considering this simple example of a actually vibrating bar which has been discretized into a set of 4 elements 1, 2, 3, 4 and the finite element discretization implies the use of these global shape functions, this is not how we actually compute though but what we do is we carry out calculation at the element level and then sum up. So what we need to do here is we need to find for example kinetic energy in terms of the displacement field as shown here and this we can add over the elements, so I can find kinetic energy contribution from this element, this will have you know contribution from only these 2 trial functions, all other trial functions are 0 here, similarly for this element there are 2 instances, for this element there are 2 instances and this is 2 instances. What we are actually doing is we are computing these 4 situations separately, so that is what we are doing, so this is what I meant when I said total kinetic energy is sum of contribution from each of the elements, so this then we can construct a Lagrangian which is again summed and if we actually then carry out the, if we implement the Euler Lagrange equation on this representation we will get basically what we got here, so we are able to achieve these calculations through this simple you know in the language of matrix algebra through these arguments. Now we will come to the question on forces, so how do we determine the nodal forces? Suppose there is a beam element which is carrying a transverse load like this, so in your finite element approximation we have approximated the field variable V of X, T in terms of U1, U2, U3, U4 and the 4 interpolation functions. Now we would like to replace this distributed load in terms of these equivalent nodal forces, so we have therefore V of X, T as I equal to 1 to 4 UI of T Phi I of X, this is how we represent the field variable and we would like to now replace this F of X, T through these equivalent nodal forces, so how I can replace therefore this situation by a combination of these two situations that's the question, so what we do is we use the principle of virtual work we define delta V X, T as a virtual displacement, we can write this in terms of virtual displacements at the nodes and the interpolation function, because this is a representation that we are talking about. Now the work done by the nodal forces on the virtual displacements must be equal to the work done by these forces on these virtual displacements, so if you now write for delta U of X, T this we can take the terms to the other side and we get this expression. Now since the delta UI of T are all arbitrary the only way this can be 0 for all T is that the term inside the brace must be independently equal to 0 and that leads us to the definition of the equivalent nodal forces, for example if the beam is carrying a UDL but as a function of time I can derive the equivalent nodal forces as P1, P2, P3, P4 as shown, similarly if you have an impulsive force acting at X equal to A I can substitute and we can get the expression for equivalent nodal forces. Now we are getting now equivalent nodal forces in the local coordinate system of the beam, how do you get in global coordinate system? It's simply the transformation, because force is a vector just like a displacement the same rule of transformation as was used for displacements also work here and this is what I get, and how do we assemble? Suppose the force vector for a built up structure, how do we assemble that? Again we follow the same logic as we use for assembling structural matrices and the global force vector in terms of element level force vector is again through a summation using these A matrices, this is as shown. So first we formulate the force in the local coordinate system and bring it to the global coordinate system, and then use the A matrices and sum over S equal to 1 to P to get the equivalent nodal forces in the global coordinate system for the built up structure. So now at this stage what we have obtained is the global structure equation of motion, the quantity that you are seeing in red here is a damping force, so I am now arbitrarily introducing this assuming that the damping mechanism is viscous and we are going to select this damping matrix in a way that facilitates the use of normal modes which are obtained by undamped free vibration analysis, therefore this C matrix can be introduced at a later stage, so I have arbitrarily introduced this matrix we will return to this issue in some greater detail later. So now the next issue that we have to consider is the boundary condition, so if you go back to this, now we have tackled a few issues, we have considered the fact that these elements carry both axial and bending energies, so we have tackled that, and the structure itself is made up of different elements which are oriented in different directions, so that we have tackled, we have first transformed all the entities in a global coordinate system and we have devised a way to assemble them, and similarly we have now tackled the issue of how to handle an external force and convert it into equivalent nodal forces in global coordinate system for the structure, at the structure level, so now the question that remains is how to deal with boundary conditions, so now, so what we do is to deal with that we partition the displacement vector into a component let us call it a UG of T, and another one let me call it a UI bar, this UG of T the G we can take it as representing the ground, the ground displacements can be specified, okay, it can be specified to be 0 or it can be specified to be a given time history, often we take that to be 0, in situations where the dynamic actions come through support motion as in case of earthquake induced vibrations UG of T will be specified time histories, so in any case they are known. The remaining nodal degrees of freedom are the structure degrees of freedom that we need to determine these are the unknown nodal displacements, now to start with let us consider that the nodal displacements are given to be 0, so the partitioning of the displacement vector will be 0 and U bar of I, so if there are Q components in UG of T this will be Q cross 1 and this will be N minus Q cross 1, now the partitioning of this displacement vector also induces a partitioning on the structural matrices, so and the forcing vector, so I will partition the mass matrix, stiffness matrix, damping matrix and the force vector in this form, this size of M00K00C00 will be Q cross Q and M2K2C2 will be N minus Q cross N minus Q, so the global of equation of motion, the global equation of motion in the partition form can be written this way. Now you see here these displacement components are specified to be 0 and they will be associated with unknown reactions and associated with unknown displacements which is UI bar there will be equivalent nodal forces which are externally applied, not all of them need to be nonzero it depends on the specific problem on hand, so there are now capital N number of equations and capital N number of unknowns, what are the capital N number of unknowns? Q reactions and N minus Q unknown nodal displacements, so we can write now the equation corresponding to the first row here, so M00 this is 0, M0I UI double dot C0I UI bar dot plus K0I UI bar is F naught of T, so this set of equations represent the equations for unknown reactions, so they can be determined only after you find this UI bar, how do I find UI bar? I write the equation from the second row, so that is this, that is MI0 into 0 plus M2 into UI bar double dot and so on and so forth the other terms on the right hand side I have the equivalent nodal forces, so these are the equations for unknown displacements, so the number of equations is capital N, Q reactions N minus Q unknown displacements. Now what we do is we avoid all this clutter of subscripts and bars and things like that, we write this equation, this set of equation which are the equations for unknown nodal displacements and velocities and accelerations simply as MU double dot plus U dot plus KU equal to F with specified initial conditions, so these are the governing equilibrium equations for the structure which we need to solve, so this has taken into account the effect of boundary conditions and all the assembling as we described, so these are the unknowns that we have to solve by solving this set of coupled second order ordinary differential equations. So now let's quickly summarize what we have been doing, we want to analyze a structure like this carrying a dynamic load like this, this structure is made up of different elements S equal to 1, I can have this is S equal to 2, 3, 4, 5, 6, so what we do for each element we first write the element level equation of motion, we know the element level mass matrix stiffness matrix and we can notionally introduce a damping matrix and we can get the force vector also for S element, then we go to the global coordinate system, okay, see the local coordinate system and global coordinate system for this element it is same, suppose if I take the global coordinate system to be this, for this element 4 this local and global coordinate system coincide not for the other elements, so you have to make a transformation, so you then get the equation of motion for elements S equal to 1, 2, P in the global coordinate system through this transformation, that means for each of the element you have to find out T1, T2, T3, TEP, they are the coordinate transformation matrix. Next for each of the element we have to now define the A matrix which connects the naming of degrees of freedom at the element level to the naming of degrees of freedom at the structure level, so that leads us to the global equation of motion after assembly of structural matrices and before imposing boundary conditions, this is the equation, this is so this is a set of capital N number of equations Q of them for reactions and N-Q for unknown displacements. We now get the equation for unknown reactions, we partition the displacement vector into Q and N-Q elements and one set of equations is equations for unknown reactions, the other set of equations is for unknown displacement, so this is the equation shown in red here is the governing equation of motion which we need to solve to obtain the solution to the problem on hand. What really happens is at the end this is a set of ordinary differential equations, T is still a continuous variable here, so these equations are known as semi-discretized equations, so we have to further develop numerical methods to discretize this time and we will get a discrete map from which we will determine all the nodal displacements, velocities and accelerations. Subsequently we have to now go back to the element level descriptions through the appropriate interpolation functions and the subsequent questions will be on how to find out displacement at points which don't coincide with the nodes, how do you find these are displacements, how do you find strains and stresses, and how do you you know compute quantities like principle stresses, principle strains, one micell stress metric etc., etc., those questions need to be still addressed.