 We are discussing questions on stability analysis of structures made of beam elements, so in the last class we derived the geometric stiffness matrix for an Euler-Bernoulli beam, and this is the total potential and we use this as a trial function, and this is the geometric stiffness matrix that we derived. We also reviewed the definition of non-linear strain displacement relations, and we derived the geometric stiffness matrix for a general solid element, and we obtained K sigma in this form, we started with an assumption that in a three dimensional solid an initial state of stress prevails in the body, and that is denoted by a sigma naught with element sigma XX naught, sigma YY naught, etc., and sigma ZX naught, and as the body strains these initial stress remain unaltered, and the work done by the these stresses on the deformation is given by this, and this gets partitioned into one set of terms which lead to the equivalent nodal forces due to the presence of the initial stress, and the other one represents the interaction between deformation and the initial stress and leading to the geometric stiffness matrix as shown here, and we using this formulation we derived the geometric stiffness matrix for an axially deforming rod, and we obtained this as a geometric stiffness matrix. Using this and combining with the geometric stiffness matrix that we derived for Euler-Bernoulli beam, we formulated the two-dimensional beam element, and that beam element to analyze frames we need to transform them into from local to global coordinate systems, and that rules for that follows what we did when we were discussing vibration analysis, and that we recapitulated, and we were asking towards the end of the class how to deal with analyzing stability of structures such as this, it could be more complicated than this in terms of loading and geometry, but right now we are considering planar frames with the discussion that we had till now, so the recipe for doing this is what we will start today's lecture with, so what we do is we assemble the reduced structure elastic stiffness matrix K and the equivalent nodal force vector F for the given structure, then we solve that problem and find out the axial thrust diagram for the structure, that is the first step, that is a simple linear static analysis, then having established the axial thrust in different elements with that as a parameter in the geometric stiffness matrix we formulate the reduced structural geometric stiffness matrix KG for each of the elements and then transform them to the global coordinate system assemble all that we will do, that leads to a KG which is a reduced geometric stiffness matrix for the structure. Then we solve the Eigenvalue problem and estimate the lowest value of lambda by solving this Eigenvalue problem, what is lambda? Suppose the loading pattern I write in a notional way as Q of XY, which represents all these concentrated forces and maybe some moments and so on and so forth, this is collectively represented by this function Q of XY, now what I do is I increase this loading pattern by a scalar parameter alpha and I will go on increasing this till the time the structure reaches a critical state, our objective is to estimate this value of alpha, so you might have designed the structure for some notional value of alpha equal to 1 for example, then Q of XY is what we are discussing here, so how much that loading need to increase before the structure reaches the verge of losing its stability is the question we are trying to answer, so that lambda that is alpha now we will estimate this lambda mean and this is a factor by which the prevailing loading pattern needs to be multiplied so that the structure reaches a critical state, so this factor is also called a load factor, the associated Eigenvector leads to a definition of the lowest buckling mode shape, so this analysis therefore has two steps, both steps are quite familiar to us, one is simply perform a static analysis, construct the axial thrust diagram, this we have done for a simple beam structures in the previous class, then using that axial thrust diagram we formulate the geometric stiffness matrix for each of the elements, transform them to global coordinate, assemble all that, after we do that we get the structure elastic stiffness matrix, the reduced structural elastic stiffness matrix and the geometric, reduced geometric structures stiffness matrix, so then we do the Eigenvalue analysis and find out this alpha is in fact lambda to get the idea about the condition under which the structure loses stability. Now I will not be discussing too many examples here because both the steps are familiar we have done something similar for both vibration in vibration analysis already and simple problems have already been solved, so I will leave this as an exercise, we will consider a frame made up of three members and which is loaded as shown here, cross sectional properties are given here, the problem is to determine the critical value of load P, then plot the axial thrust diagram at the critical loading condition and plot the corresponding buckling mode shape, so you could use three elements, one element for this, another element for this and another element for this, if you go back in your lecture notes, one of the earlier lectures where we studied the vibration analysis we have analyzed the same structure, so the stiffness matrix, elastic stiffness matrix is already available, you have to follow the same steps but assemble the geometric stiffness matrix and perform the Eigenvalue analysis, now in the discussion so far we have discussed about determining the load factor that helps us to determine the possible existence of a neighboring equilibrium state in the structure, but we could also do a response analysis by including geometric stiffness even when the structure is not on the verge of buckling, if the structure when the loading may not have reached a critical state but still I may like to include effect of geometric stiffness in characterizing the response, so to be that we have done in the study of beam columns, we got the stability functions and how the response amplifies all that we have discussed, so to do this now total potential we'll write this is the work done by the loads, and again use the trial function such a like this and construct the function U and write the equilibrium equation, here I will get the equilibrium equation to be K-PJ into U is equal to Q, and U is obtained by inverting this relation after imposing suitable boundary conditions etc, similar exercise can also be done for this type of problem also, even when P is not close to a critical state we can find out the influence of geometric stiffness on the response, as an exercise what we could do is we can revisit this beam column problem, we have considered a simply supported beam with a concentrated load Q and actual load P, and when the load is applied at the mid span and if you are focusing on the mid span displacement we have found out the mid span displacement to be the displacement with P equal to 0 multiplied by this chi of U which is the stability function, U is lambda L by 2 lambda itself is square root P by EI that we have done earlier, so we have these three stability functions chi, epsilon and xi, chi characterizes displacement, epsilon characterizes the rotation at left end and M this xi characterizes the amplification with the bending moment, now what we could do is you can carry out an analysis using finite element method and compare the response obtained using this exact calculation, analytical calculation with the finite element solution, we see the influence of number of, you know choice of number of elements and such issues can be examined using that, this is an exact solution that is available, so this forms a firm benchmark against which you can check your answer, so this I leave it as an another exercise. Now that is effect of geometric stiffness on static behavior of the structure, now we could also do a similar analysis using include geometric stiffness as well in dynamic analysis, so this we will consider slightly later we will study the effect of axial loads on natural frequencies and force response of structures but as a prelude in this context it is worth mentioning that this is also possible, so this generalization of beam column type of analysis for dynamic loads, here the axial load modifies the natural frequencies and mode shapes and naturally the force response also would get modified, now we have considered so far the coupling between the axial load and bending, so we can also ask the question if there is coupling between axial load and twisting, okay, so see the issue is we are using non-linear strain displacement relations, so a given strain component is related to several gradients of UVW, so a coupling has to be expected, I mean that should not come as a surprise, so to be able to do that what we do is we recall we have discussed that problem of torsion of shafts in lecture 11, so from there we pick on, we are considering now XZ plane and we assume displacement V-Z theta X and WY theta X and if we use now the strain of a fiber along longitudinal direction based on non-linear strain displacement relation we get this, now focusing on geometric effects we have to include only the non-linear terms, the linear terms contribute to equivalent nodal forces arising out of pre-initial existing stresses that is not we are interested in, so we compute the strain component, this G is geometric and this is what we are interested in and using the assumed displacement field we get this to be given by this, now the only pre-stress that, you know initial stress that exists is sigma XX naught due to axial load, so the question we are asking is due to this axial load we normally expect the beam to deform along X axis, but for some value of P it can twist about X axis, that is what the value of P which induces a twisting in the beam is what we are computing, in the shaft that we are computing, earlier we have computed the value of P which induces bending, right, now the twisting, so that also can be done, so this stress sigma XX naught acts on this longitudinal strain leading to this term and now we can, where sigma XX naught is, we write P by A for sigma XX naught and I P is this area moment of inertia, now we can use finite element method, the first order derivative is there and two-noded element means theta X should be the degree of freedom, use linear interpolation functions, we have done this several times, so the geometric stiffness arising out of this interaction is given by this, okay, so now what we have done is interaction between axial effect and planar bending and coupling between axial load and twisting, so with this we seem to be ready to deal with problem of a three-dimensional beam element. Now if you recall what we did for assembling elastic stiffness matrix, assuming that the beam cross section has some kind of bi-directional symmetry, there is not an symmetric cross section because of which there will be coupling between bending and twisting, etc., those conditions we are not considered, nor are we going to consider that here also, assuming that kind of symmetry exists in formulating elastic stiffness matrix it was simply a matter of assembling the contributions from the two bending actions, twisting action and axial deformation action, there was no new phenomena that came in because we considered a three-dimensional element, but now we need to ask this question when it comes to discussing the effect of geometric stiffness or of characterizing the geometric stiffness, is it simply a question of assembling and just instead of dealing with a lower order matrix we deal with a higher order matrix or are there any new phenomena, okay. So we can start with a very simple model, we can include effects of coupling between axial forces with bending along the two axis and torsion about the longitudinal axis, so in a three-dimensional beam due to axial load the beam can buckle in two directions, for example if beam cross section is this, one value of P will induce buckling along this direction, another value of P will induce buckling in this direction, and similarly some value of P will induce buckling in twisting mode, so all these three coupling effects we have the coupling between axial load and bending here, bending here and twisting we have included, so equipped with that now we can simply assemble the matrix, including these effects, so this is if you name the degree of freedom as shown here two nodes and this is a nomenclature of degrees of freedom, the field variables are U, V, W and theta, if you do that if you recall this we already seen these for example 1, 1 element will be due to axial deformation, similarly 7, 1 will be due to, see for example 7 and 1, 1 and 7 are axial deformation, so etc., so this will get filled up, and if you write the elastic stiffness matrix we have seen this matrix earlier, this is the elastic stiffness matrix, so AE by L, minus AE by L, minus AE by L, AE by L are here, and similarly you can identify the other elements for bending in XY plane, XZ plane and also twisting about theta X, this is the story as far as elastic stiffness is concerned, the similar story repeats now for the geometric stiffness matrix, so there is again axial effects bending about two planes under twisting effect, these effects are included, so this seems to start with a reasonable model to work with, but one should understand that this is not enough, even for a simple model, the initial state of stress in the beam actually could include presence of longitudinal force, torque along longitudinal axis and two bending moments, so when we are talking about initial stress why think of only axial loads, there will be initial stress is 3 dimensional, 3 dimensional beam element has several stress resultants and the state of stress is 3 dimensional, so there will be a longitudinal force and two bending moments and twist along the longitudinal axis, these four stress resultants will be present, so in the initial state of stress all four will be there, so the question arises how they will interact, why only axial loads should interact with bending, why not bending interact with torque or torque interacts with another plane or deformation so on and so forth, so a general formulation within the framework of Euler Bernoulli theory need to include interaction of torque and bending moment in addition to axial force and spatial variation of torque and bending moments within an element, so it is reasonable to assume axial deformation is constant over an element, but bending moment and twisting moment need not be constant, so because in a pre stress when you draw your bending moment shear force and axial thrust and torque diagrams they need not be constant, they can vary within an element, so of course other complications are asymmetric cross sections, this we are not considering amnesotropy and such things, so we will restrict to the sections which are bi-symmetric and try to include these effects, the fact that such effect should be present again I want to emphasize should not come as a surprise because if you see the structure of the expression for U in which the work done by the initial stress on the deformation is characterized there is quite a wide range of coupling terms present here, so for a three-dimensional beam element moment based on Euler Bernoulli beam hypothesis we formulate the expression for UVW, we can fill up all these terms, right, you assume certain deformation field in Euler Bernoulli beam formulation, you fill it up here, so you will see which interaction is present and things like that, so if we do that which I am not going to get into all the details, we will get a more populated geometric stiffness matrix, this for example has been derived or discussed in this book by McGuire, Gallagher and Zeeman, I have taken the information from this book it is given here, okay, so again let me emphasize when we come to three-dimensional beam elements, when we formulate geometric stiffness matrix there are newer phenomena which we need to be concerned about because of complicated couplings that can exist between bending, twisting and actual, you know, deformation. Now we have in our study on dynamics we have studied plane stress element, plate bending elements, triangular element, quadrilateral element, isoparametric formulation, etc., etc., solid elements so on and so forth, all those elements now we can revisit and derive the geometric stiffness matrix, so we have to consider for a given problem the appropriate form of strain displacement relations which need to be nonlinear now, and identify all the coupling terms and I derive the K sigma matrix, that can be done. Now what I will do is, I will not do all the calculation for all the elements that we have done, that exercise is something that can be done, but however I would like to consider this problem of geometric stiffness matrix for a thin plate bending element following Kirchhoff-Lauw hypothesis. So this plate has the initial stresses which are planar that is the membrane forces NY, NX, Y, NX. Now the membrane forces are results of the stresses sigma XX, sigma YY and sigma XY, which we typically carry out using plane stress analysis, and the stress resultants NX, NY and NXY are obtained by integrating these stresses across the depth as shown here. Now we will assume that in a plate problem there is a prevailing membrane stresses initially, and that structure bends. So when the straining takes place that is due to bending, the initial stress that exist, that is the membrane initial stress that exist will not change, that is assumption that we are making. So the idea is to examine the work done by constant membrane forces on the displacements associated with transverse deflections. So membrane strains associated with small rotations, you know these are given by this, okay, this is extract from the appropriate strain displacement relations. I am retaining only the relevant nonlinear terms, the linear terms lead to equivalent nodal forces associated with initial stress that I am not discussing. Now we can assume that NX, NY, NXY are independent of W, that means as W the transverse displacement occurs these things, these NX, NY, NXY do not change. Now the work done therefore is given by this. So now for epsilon XX, epsilon YY, and epsilon XY I use these nonlinear strain displacement relation terms from that, and I get this equation, okay. Now again here W is the field variable and using the matrix of shape function N and nodal coordinates UE I will derive this, okay. This leads to the G matrix that we have seen earlier, and W is written as G into UE. Now I compute K, the elements of K sigma matrix using this, this relation we have already derived and now the integration goes over the area. I am not giving all the details, just giving you the main you know milestones in the derivation. Now in a isoparametric formulation what we do, we need to you know introduce the coordinate XI and eta, and this evaluation of this integral has to be based on this transformation, and that leads to this Jacobian and a different matrix on a GI which is now relates DOW by DOXI and DOW by DO eta to the nodal coordinates. This is the GI matrix. Now the K sigma is given by this, where the integration is done over XI and eta which run from minus 1 to plus 1. Now if membrane forces are taken to varying space within an element that can be incorporated here in evaluating this, but then while choosing order of integration you should have an idea what is the order of variation of this NX, NY, NX, Y. Typically you can assume linear or quadratic variation and then decide upon the order of integration. K sigma do not depend on elastic properties of the plate, see if you see here there is no Poisson's ratio and NX modulus here, but of course the P stress, the initial stress that exists will be function of E and sigma, E and mu, but this elements of this itself are not directly dependent on elastic constants. Therefore this formulation can therefore be used if material is unisotropic also, because this doesn't depend on material being isotropic. Now we can generalize this formulation to thick plates, as an illustration we can consider three dimensional Timoshenko beam, so the game is now to write the correct expression for U sigma. So in Timoshenko beam theory you have a postulation for the displacement field, based on which you can compute the strains and using the general expression for U sigma and retaining the appropriate terms you can always derive the appropriate functional. Now I have just listed out the functional here for Timoshenko, three dimensional Timoshenko beam in terms of N bar the axial force, MX bar, MY bar and MZ bar which are the stress resultants in the initial state that prevails in the three dimensional beam element, and this is a work done by deformation of the beam on this initially prevailing stress field. So this we can use now a suitable interpolation scheme and evaluate the K sigma matrix. So this again if you see look at the functional the highest derivative is 1 and the field variables are V, W and theta Z, so you can assume linear interpolation functions for these three run through the whole formulation you will get K sigma. So all these are following, setting up this expression for U sigma all other steps are we are already familiar with, it's only matter of doing them correctly that will take you to the final answer. There is a reference that you may find useful, this is a handbook of mechanical stability in engineering and this discussion on three dimensional beam is available in volume 2, stability of elastically deformable mechanical systems. Now we will now move on to discussing few conceptual issues associated with stability of structures, and I want to now discuss what are known as imperfection sensitive structures, what are they? To be able to do that we will not get into elastic systems but this idea can be explained by considering rigid link which is rigid bar which is loaded by supported through springs as shown here. There are three arrangements, in the first arrangement the bar is hinged here and supported on a rotary spring, in all the arrangement the bar is hinged at A, and in the second arrangement the bar is supported at the point B through this spring, whereas in third case the bar is supported through an inclined spring along an inclined axis as shown here. Now I will call this as configuration 1, 2 and 3. Now we can perform a thought experiment now, the thought experiment is such that let us consider we are selecting K 1, K 2, K 3 for these three problems in such a way that the critical load for all these three systems is identical, they are all equal, so PCR 1 equal to PCR 2 equal to PCR 3, so this K 1, K 2, K 3 are adjusted in such a way that all these three structures should lose their stability at same value of P. Now design 1, 2 and 3 in such a way that PC 1 is equal to PC 2 equal to PC 3, that is select the spring constant K 1, K 2 and K 3 such that these three critical loads are equal. We can verify that for the first system the critical load is given by K 1 by L, second one it is K 2L, third one it is K 3L cos square alpha, where alpha is this angle, so you should notice that the units of K 1 and will be different from those of K 2 and K 3. So what we are doing is we are selecting now K 1, K 2, K 3 such that K 1 by L equal to K 2L and that is equal to K 3L cos square alpha. Now let's assume that this is a thought experiment, we are fabricating these three systems carefully and performing this analysis, load deflection analysis. No matter what care we take in fabrication and no matter what care we take in applying this load in an ideal manner, it will always be observed that PC 1 will be greater than PC 2 and that will be greater than PC 3. Although as per the design all three loads must be equal, not only that if you were to repeat this experiment several number of times the scatter that you get in this third case will be more than the scatter in the second result and that scatter will be more than the scatter in the first result. So this apparent you know looks like a paradox, so what we are doing is if we were to fabricate the three models and perform experiments it would be observed that PC 1 will be greater than PC 2 and that will be greater than PC 3, no matter how carefully the three models are fabricated and no matter how diligently the experiments are performed, why it happens? It also would be observed that the scattering the experimentally determined value of critical loads would be such that the scatter in the third system will be more than the scattered in the second system and that would be more than scatter in the first system. Why this happens is a question that we wish to now answer, in fact this analogy is not a very theoretical one, we can show that a beam structure behaves like this, a plate structure behaves like this, a shell structure behaves like this, so that would mean if you select a beam, plate and a shell structure to have a same critical load you want, if you perform an experiment the observed critical load will be you know for a beam will be greater than that of a plate and for a plate it will be greater than that of a shell, why it happens is a question that we are trying to answer. Now let's analyze this problem, so let's assume an ideal situation where every load is applied perfectly the structure is supported perfectly and theta is a degree of freedom, if you assume that this neighboring equilibrium state exists then V is given by half K theta square minus PL 1 minus cos theta, now equilibrium dou V by dou theta must be equal to 0 and that gives the load deflection diagram which is P is equal to K theta by L sine theta, I am not assuming theta to be small, for stability for every value of pair of P and theta here if this K minus PL cos theta is greater than 0 that point would be stable on the load deflection diagram otherwise it is unstable. Now this is an ideal situation, we will introduce a notional imperfection in the system, let's assume that there is some kind of free play in the support here and the spring will come into action only after theta crosses an angle epsilon which is small, so the V in this case will be half K theta minus epsilon whole square and this, again for equilibrium I get this equation and for stability I get this equation, the stability condition is independent of epsilon but the load deflection diagram depends on epsilon, now if I plot this, this is how it will look like, suppose epsilon equal to 0 this is one Y axis I have the load and on X axis I have theta, so as the load increases the load deflection path rises along this curve and once the P critical value is reached the theta equal to 0 solution becomes unstable, red lines are unstable and black lines are stable, where in the stable branch bifurcates into two other branches like a bowl as shown here, okay. Now in presence of an imperfection, suppose epsilon is pi by 50, a small imperfection the load deflection path has a rising part, it you know climbs up like this and it goes along this, and there is another branch where initially it will be unstable but afterwards it is against stable, okay, so that would mean at this value of theta this is stable branch. Now we can make some observation that the load displacement curve for imperfect system rises along a stable path without encountering any instability, right, it won't encounter any instability, then the deflection grows more rapidly as the load approaches the critical value corresponding to a perfect system, see if a linear analysis is done if you recall our answer will be like this, right, this is the loading path, this is the horizontal, the nonlinear analysis takes it to a more realistic formulation. As you come in an imperfect system, as you come closer, even perfect system, as you come close to the critical load of elastic system, sorry linear system, the rate of change of load is rapid, I mean rate of change of deformation is rapid with increase in load, that is how the kind of remnant of instability is manifested here. This is an example of a symmetric stable bifurcation, beam elements display this type of behavior, later on depending on the time available we may discuss problem of an elastica and we will be able to show that something similar is observed in a beam. Now let's consider the second configuration where the bar is supported on a spring as shown here, again the potential is given by half K L sine theta whole square minus P L 1 minus cos theta, for equilibrium this is the condition, this is a in fact is a load deflection diagram, there are two equilibrium paths sine theta equal to 0 and this equal to 0, now for stability this condition must be satisfied, okay, this is straight forward derivation. Now again let us assume that there is a free plane and the spring remains, let the spring remain unstrained till theta crosses epsilon, then V is given by this and the condition for equilibrium is this, this is a load deflection path and this is a condition for stability. So we can write a small program, plot this and at every value of P, theta you substitute into this and see whether the point is stable or not, that is what I have done. Now if you plot this something interesting happens, the load for perfect system the load deflection theta equal to 0 is the stable equilibrium path, it rises like this and then it bifurcates into two unstable paths, this is unstable path, okay. Now in presence of imperfection how does its system behave? The load deflection diagram the load deflection path increases but after it reaches a value it encounters a point of instability, so beyond this point it cannot carry further load because that state becomes unstable, okay, so the load carrying capacity is now inhibited by presence of imperfection, whereas here that was not so, it continues to carry the load but the rate of change of deformation changes as you come close to the critical point, whereas here it simply ceases to be stable and if there is greater amount of imperfection instead of pi by 50, if it is pi by 20, this is the critical load now load carrying capacity is this, it reduces, okay. So the observation here is that the load displacement curve for imperfect system is not always stable, instead it rises along a stable branch and exhibits a limit point beyond which the path drops and becomes unstable, greater the imperfection lesser is the ability of the structure to carry the load, this is an example of a symmetric unstable bifurcation, so the bifurcation is symmetric, the two paths that emerge after bifurcation are symmetric with respect to this line and a stable path becomes two unstable path, so it is symmetric unstable bifurcation, plate structure display this type of behavior. Now let's come to the asymmetric structure, suppose this is the bar which is supported on spring like this at some time, at some snapshot it would have deformed like this and we have to set up the expression for the total potential, so we have to find out the amount of stretch in the spring and the work done by gravity, so I mean the movement through this distance, so if initial length is OB is S and OB bar is S bar, the strain energy stored in the spring is half K S bar minus S whole square and this is minus PL cos theta, now using geometry I can find an expression for S and S bar in terms of this theta, theta is a degree of freedom, so if you carry out a simple calculation we will find out S and S bar and if I substitute into this, this is the expression for V, again equilibrium is dou V by dou theta equal to 0 and this is a load deflection path, now more involved and for stability I have to differentiate this further and see the sign of the term. If you perform this calculation for a, when there is no imperfection the load deflection path rises along a stable path and it bifurcates with one stable path and another unstable path, so the unstable path is of greater interest to us because that is what we are worried about, now in presence of imperfections what happens is the load increases and the load displacement path increases and again encounters a limit point afterwards it becomes unstable, there is of course a stable path along with the load can go on increasing but that is not what governs our design, this branch governs the design, so here again the stable equilibrium path is asymmetric with respect to the undeflected configuration state with respect to which the bifurcation takes place, okay, and that is actually the structure itself is unsymmetric, so you can expect something like that to happen, the stable load displacement curve for imperfect system is not always rising and stable, instead it rises along a stable branch and exhibits a limit point beyond which the path drops and becomes unstable, greater the imperfection again lesser is the ability of the structure to carry the load, this is an example of an asymmetric unstable bifurcation, shells and some frame structure display this type of behavior, so in summary we have now three situations where presence of imperfection lowers the load carrying capacity of these two type of structures, so we call them as imperfection sensitive structures, so the maximum load carrying capacity for these two systems can be shown as a function of epsilon and it drops off like this, greater the imperfection lesser is the ability to carry the load and here again there are one unstable branch and a stable branch we can plot this and you can actually compute these graphs exactly, although I have shown it as schematic here, so this clarifies the meaning of word imperfection sensitive structure. Now these three examples that we considered actually they together illustrate the possibility of diverse type of behavior which could be displayed by different structures, if a structure differs the function U, the function U changes and U as a function of generalized coordinates is different for different types of structures, so these three examples display the dramatic changes that can happen in the behavior of the structure due to the dependence of U on the generalized coordinates, the function U of theta encapsulates different types of structures, U theta can be expanded in Taylor series about the equilibrium point and different bifurcation patterns can be discerned by examining the signs of various terms. This approach forms the basis of Quiter's theory for stability analysis, I am not going to get into the details, the book by Thomson and Hunt has good description of this and monograph on Quiter's work itself is available that can also be referred to, here the non-linear analysis helps us to examine the likely importance of geometric imperfections in a given context, so it is a interaction between non-linearity and the geometric imperfections that explains the paradox that I mentioned, it is because that imperfection, no structure can be fabricated with, in an ideal manner the imperfections are ineligible, a beam structure is, performance of a beam structure is relatively insensitive to that imperfection, whereas these two structures due to presence of imperfection the load carrying capacity is lowered and the amount by which the load carrying capacity gets lowered depends on the actual, there actually one can work out that detail as indicated here, so the paradoxical observation that we started with, the explanation lies in examining the interaction between imperfections and large deformations of the structure. Now we will now consider a few more situations where questions on stability of structure arose out some surprises, so we will start with discussion on beam on elastic foundation, so let's consider a simply supported actually loaded beam on Winkler's foundation, the basic assumptions in Winkler's foundation is the reactions are proportional to the deflection at any cross section, and the spring works both in tension and compression, and actually the foundation behaves like a series of springs without mutual interaction, that is a basic assumption of Winkler's foundation. So the expression for U in this case has these two terms which you already encountered, this is a new term, half KY square DX and this is the work done by external loads which we need to include to completely formulate the problem. Suppose if you assume P of X to be 0, we get this field equation and these are admissible boundary condition, you can put it in the Hamilton's principle and derive this and derive the appropriate boundary conditions. Now as I said let's consider a simply supported beam on elastic foundation, which carries axial load, so the governing equation will be of this form. Now these are the boundary conditions at X equal to 0, displacement and bending moment are 0, and X equal to L, displacement and bending moment are 0, so for purpose of discussion we will assume that the solution can be expanded in sine terms, and I substitute into this equation and I get this characteristic equation and the P critical, that means the value of P for which a non-zero solution is possible is given by this, okay. This follows the way we analyzed Euler buckling load for a simply supported beam, if we got this P critical, if K is 0 we got this result earlier, okay, M square pi square by L square into EI. Now you should notice now the dependence of P critical on M, and M is appearing in the denominator here, so this creates some surprises, that would mean if M equal to 1 need not lead to the lowest critical load, okay, depends on values of K and EI, so what that means? Suppose if I now want to find out the value of M at which P critical is smallest, so I will express P critical as a function of M, and suppose if I now use the criteria D by D M, P critical M equal to 0, then I get an equation, and from this I get the relation between M and all these parameters, if M satisfies this equation that is the value of M for which P critical will be the lowest. So if I substitute that optimal value of M into the P critical value I get this as 2 into square root K EI, this value is independent of M. Now that means it is possible that P2 is less than P1, suppose P1 is a load at which we get a half sine wave, and P2 is a load at which we get a full sine wave, it is not necessary that the ordering should be such that P2 is greater than P1, P2 can be less than P1, so this can be the critical load depending on a problem that you think of. So that is plotted here, if I plot P critical versus L by R whole square on a log log scale, what we see is for different values of M, 1, 2, 3, 4, we see that all these squares you know don't fall below this value of 2 into square root K EI, that is this value that we derived, and if your L by R ratio is somewhere here then the lowest critical value corresponds to M equal to 1, but on the other hand if L by R ratio is somewhere here you see that the black curve M equal to 2 will be the lowest critical value, the next one will be M equal to 2 and 3 seem to have the same value, okay, so this is something that you know is special to beams on elastic foundation. So now as an exercise what we could do is we can consider an Euler Bernoulli beam element resting on a Winkler's foundation, and you can develop the elastic and geometric stiffness matrices for the element, now using the element thus developed determine the critical load for a simply supported beam and illustrate the result shown in the figure on the previous slide, that means what I want is you consider a single span simply supported beam on elastic foundation formulate the elastic stiffness and geometric stiffness matrices, assemble them by you choosing you know desired number of elements and try to see whether you get this type of behavior, whether the finite element model captures this and how good is the finite element model because these are exact results, I leave that as an exercise. Now we can consider certain other problems for example buckling of a ring, suppose I have a ring with radius R, so AB is a deformed geometry of the ring, O is the center, so OM will be R and OM1 is the deformed configuration, so MN is a section of length DS in the under deformed configuration, so MN and in the deformed configuration becomes M1N1, now DS is nothing but RD theta and therefore I get D theta by DS is 1 by R, now upon deformation the section MN becomes M1N1, now let rho be the radius of the deformed configuration, so by following this argument 1 by rho will be D theta plus delta D theta divided by DS plus delta DS, so DW by DS consequently the angle between tangent to a center at M1 and the perpendicular to the radius MO that is slope basically, if DW by DS is angle at M and at N I will have DW by DS plus the gradient into that D square W DS square into DS, so therefore based on that I get delta D theta as this, so similarly if you consider length M1N1 it is delta DS which is given by this and we can show that this is, you can work through this and you can get this equation, so you can thus get the expression for 1 by rho which is the radius of curvature for in the deformed configuration, assuming these terms here and simplifying and retaining only first order terms after doing as expansion I will skip this discussion you can verify that it is correct, we get this equation, okay. Now if the beam ring is now loaded by a uniform force field as shown here, can it assume a neighboring equilibrium state, for example as shown in this elliptical curve, this is the ring in the undeformed configuration, Q is the load which is axisymmetric constant external compressive load, so how will the ring deform? Now can we find out critical value of Q at which a neighboring equilibrium becomes possible, so we can consider the equilibrium equation and write the expression for bending moment we can do that and I get the governing equation to be given by this, there QR is a hoop force, you can do a simple calculation and get this equation, this will be the governing equation and if I now assume W of theta is A cos K theta plus B sin K theta plus this particular integral, now the DW by D theta will be minus A K sin K theta plus B K cos K theta. Now the boundary conditions that we use is for the shape that we are looking for, this shape DW by D theta is 0 at theta equal to 0 and pi by 2 by simple symmetric argument, if I do that this B becomes 0 and the second condition leads to this characteristic equation and based on that I get the value of QC to be 3EI by RQ, okay, so this is a simple of how we can compute critical load for a ring, now this is an exact solution, so what we can do now is you can consider a ring and assume Euler Bernoulli beam elements I can approximate the circle by straight beam elements, maybe use maybe 12 or 18 number of elements and compute formula, we have now the geometric stiffness matrix and elastic stiffness matrix, so you need to transform them correctly, assemble them and if you carry out the analysis for critical condition you will get an approximation and the exercise is to analyze the approximate solution we saw with this exact solution, so that is left as an exercise. So with this we will conclude this lecture, in the next lecture we will consider some questions about arch buckling and then move on with other topics in stability analysis.