USFDC Home  USF Electronic Theses and Dissertations   RSS 
Material Information
Subjects
Notes
Record Information

Full Text 
xml version 1.0 encoding UTF8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd leader nam 22 Ka 4500 controlfield tag 007 crbnuuuuuu 008 s2010 flu s 000 0 eng d datafield ind1 8 ind2 024 subfield code a E14SFE0004558 035 (OCoLC) 040 FHM c FHM 049 FHMM 090 XX9999 (Online) 1 100 Shitta, Abiola. 0 245 Modeling swelling instabilities in surface confined hydrogels h [electronic resource] / by Abiola Shitta. 260 [Tampa, Fla] : b University of South Florida, 2010. 500 Title from PDF of title page. Document formatted into pages; contains X pages. 502 Thesis (MCHE)University of South Florida, 2010. 504 Includes bibliographical references. 516 Text (Electronic thesis) in PDF format. 538 Mode of access: World Wide Web. System requirements: World Wide Web browser and PDF reader. 3 520 ABSTRACT: The buckling of a material subject to stress is a very common phenomenon observed in mechanics. However, the observed buckling of a surface confined hydrogel due to swelling is a unique manifestation of the buckling problem. The reason for buckling is the same in all cases; there is a certain magnitude of force that once exceeded, causes the material to deform itself into a buckling mode. Exactly what that buckling mode is as well as how much force is necessary to cause buckling depends on the material properties. Taking both a finite difference and analytical approach to the problem, it is desired to obtain relationships between the material properties and the predicted buckling modes. These relationships will make it possible for a hydrogel to be designed so that the predicted amount of buckling will occur. 590 Advisor: Ryan Toomey, Ph.D. 653 Buckling Polymer Elastic Stability Finite Difference Approximation Biharmonic Equation 690 Dissertations, Academic z USF x Chemical & Biomedical Engineering Masters. 773 t USF Electronic Theses and Dissertations. 4 856 u http://digital.lib.usf.edu/?e14.4558 PAGE 1 Modeling Swelling Instabilities in Surface Confined Hydrogels b y Abiola Shitta A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Chemical Engineering Department of Chemical and Biomedical Engineering College of Engineering University of South Florida Major Professor: Ryan Toomey, Ph.D. Babu Joseph, Ph.D. Daniel Simkins, Ph.D. Date of Approval: July 1 2010 Keywords: Buckling, Polymer Elastic Stability Finite Difference Approximation, Biharmonic Equation Copyright 2010, Abiola Shitta PAGE 2 Dedication I would like to dedicate this work to my father and mother to whom I am and always will be eternally grateful They took the time to make sure that I had all the tools necess ary to become the best that I could be. Had it not been for this, my life would be very different and I probably never would have lived up to my greatest potential . PAGE 3 Acknowledgements I would to thank Dr Toomey for allowing me to join his research group and also for allowing me to pursue a topic that was of great interest to me. I would also like to thank Samuel DuPont and Ryan Cates for their insight through discussion as well as through the data that they collected in their own work. I am greatly thankful for the help and assistance offered by Ali Gardezi and Mario Juha. Ali took the time to impart his depth of knowledge of MATLAB on me. Without his help, I would not have had the confidence to pursue this topic. Mario provided me with a lot of the background of the theories governing Continuum Mechanics and Elastic Theory. This understanding was paramount and without it I would not have been able to finish this thesis in a timely manner. Finally, I am thankful for t he direct and indirec t support of all the other members in the research group. PAGE 4 i Table of Content s List of Figures ................................................................................................................... iii Abstract ............................................................................................................................ iv Background ....................................................................................................................... 1 Continuum Mechanics ........................................................................................... 2 Linear Elastic Theory ............................................................................................. 4 Structural or Elastic Stability Theory ...................................................................... 4 Biharmonic Equations and the Airy Stress Function ............................................. 5 Derivation of Airy Stress Function ......................................................................... 6 Finite Difference Approximation ........................................................................................ 8 Basic Idea and Setup ............................................................................................ 9 Stability of Solution .............................................................................................. 10 Boundary Conditions ........................................................................................... 11 Matrix Inversion as a Means of Solution.............................................................. 13 Analytical Solution ........................................................................................................... 15 Separation of Variables ....................................................................................... 15 Using Boundary Conditions to Solve for Constants ............................................. 19 Determination of Critical Buckling Values ....................................................................... 23 Failure of Plane Strain Assumption ..................................................................... 23 Determinant Method for Critical Buckling Values ................................................ 24 Energy Method .................................................................................................... 27 Bending and Wasted Modes of Buckling ........................................................................ 30 Equations Involved with Different Buckling Modes .............................................. 30 Possible Reasons for Different Buckling Modes .................................................. 32 Results ............................................................................................................................ 34 Finite Difference Solution .................................................................................... 34 Analytical Solution ............................................................................................... 37 Determinants ........................................................................................... 38 Critical Values .......................................................................................... 41 Conclusion ...................................................................................................................... 46 References ...................................................................................................................... 48 PAGE 5 ii Appendices ..................................................................................................................... 50 Appendix A: Nomenclature .................................................................................. 5 1 Appendix B: Derivations of Finite Difference Equations ...................................... 52 Appendix C: Derivation of Biharmonic Equation in Terms of Displacement ........ 54 Appendix D: Proof of Solution to Biharmonic Equation ....................................... 56 Appendix E: Derivation of the Strain Energy Equation ........................................ 58 Appendix F: Determinant of a 4x4 Matrix ............................................................ 60 Appendix G: MATLAB Code ................................................................................ 62 PAGE 6 iii List of Figures Figure 1: Setup of Buckling Problem ................................................................................. 2 Figure 2: Boundaries and Grid Definition of Finite Difference Method .............................. 9 Figure 3: Boundary Conditions at Specified Boundaries ................................................. 13 Figure 4: Determinant Behavior w.r.t. Nx ......................................................................... 25 Figure 5: Plot of Solution of Critical Values of Nx ............................................................ 26 Figure 6: Work and Strain Energy Plots Used for Critical Values of Nx ........................... 29 Figure 7: a) Wasted and b) Bending Modes of Buckling ................................................. 30 Figure 8: Buckling Transition from One Mode to Another ............................................... 32 Figure 9: Spy of the Coefficient Matrix of Finite Difference Problem ............................... 35 Figure 10: Solution Instabilities with Different Numbers of Grid Points on the y Axis. .... 36 Figure 11: Solution Instabilities with Different Numbers of Grid Points on the x Axis .... 37 Figure 12: An Analytical Solution with Five Half Waves .................................................. 38 Figure 13: Determinant Behavior w.r.t. Nx at Different Values of m ................................ 39 Figure 14: Determinant Behavior w.r.t. Flexural Rigidity at Different Values of m .......... 40 Figure 15: Determinant Behavior w.r.t. Length and Width Ratio ..................................... 41 Figure 16: Critical Values of Flexural Rigidity at Different Values of m and Nx ............... 42 Figure 17: Critical Values of Nx at Different Values of m and Flexural Rigidity ............... 43 Figure 18: Critical Values of Nx at Different Values of m and Length to Width Ratio ...... 44 Figure 19: Critical Values of Length to Width Ratio at Different Values of m and Nx ...... 45 PAGE 7 iv Modeling Swelling Instabilities in Surface Confined Hydrogels Abiola Shitta Abstract The buckling of a material subject to stress is a very common phenom enon observed in mechanics. However, th e observed buckling of a surface confined hydrogel due to swelling is a unique manifestation of the buckling problem The reason for buckling is the same in al l cases ; there is a certain magnitude of force that once exceeded, causes the material to deform itself into a buckling mode. Exactly what that buckling mode is as well as how much force is necessary to cause buckling depends on the material properties Taking both a finite difference and analytical approach to the problem, i t is desired to obtain relationships between the material properties and the predicted buckling modes These relationships will make it possible for a hydrogel to be designed s o that the predicted amount of buckling will occur. PAGE 8 1 Background Responsive pol ymers a nd their applications have become a very active area of research. The applications for these polymers range from biomi metic actuators to biocatalysts and the list continues to grow [1]. R esponsive surfaces are another application of these materials where the fact that the polymer is attached to a surface can give rise to unique behavior [2, 3]. A responsive polymer hydrogel confined to a surface can be made to swell with a variety of different stimuli. This swelling can, under certain circumstances, cause the hydrogel to buckle. If th e gel is free to swell as it pleases or is not constrained, there is no observed buckling phenomenon. However confining the hydrogel to a surface causes stress to build up within the gel as it swells. This build up in stress is caused by the opposition to its free swelling due to the surface confinement. This opposing force manifests itself as a compression of the gel. If the magnitude of the compression is sufficient, buckling of the hydrogel will occur. Fig. 1 shows an illustration of this system with a uniaxial compression. This is the framework that will be used to model buckling in the hydrogel. PAGE 9 2 Confining Surface Hydrogel x y z Compression Compression Figure 1 : Setup of Buckling Problem The stress at which buckling occurs depends on a variety of factors; the dimensions of the m aterial, the way it is supported and the properties out of which it is made [4]. In the case of a swelling hydrogel, it could also be dependent on how much the hydrogel is al lowed to or capable of swelling since this could very well determine the magnitude of c ompression that it undergoes. It is the goal of this paper to find out how these factors play a role in the buckling of a hydrogel and potentially what relationships can be determined from these phenomena. Continuum Mechanics Continuum mechanics make it possible for a body being subject to stress to be analyzed with minimal complexity A continuum body is described as a body in which the molecular interactions are ignored. If yo u were to continually break the body up into smaller and smaller pieces there will be a point at which the interactions between those pieces will be negligible. Those pieces w ill also eventually become approximate points which cover the entire volume of the continuum body. This idea makes it possible to solely analyze the effect of applying a load to a material without the added complexity of studying the intermolecular forces. PAGE 10 3 Within continuum mechanics there are principles that apply to all materials as well as the constitutive equations that define specific mechanical behavior. T he general principles refer to things like the conservation of mass, energy and momentum. To our knowledge, all known materials obey these general principles; so the equations that they generate set up the foundation of continuum mechanics. The constitutiv e equations apply to materials that are linear elastic isotropic, linear elastic anisotropic, nonlinear elastic isotropic and so on and so forth. These descriptions tell exactly what modifications have to be made so that the general principles still apply and accommodate the behavior of the material. The difficulty with problems in continuum mechanics, as with problems in any other field, lies primarily in how much information you have about the problem at hand. For example, let us say that you know the ini tial and final conformation of a material and you want to find the equation s that describe what stress es are required to make that so. This problem would require that you have an equation that describes the displacement of the points in the body from the i nitial to the final state. This equation is difficult to come by but it could probably be determined by marking the material initially and observing how the same points move after loading to arrive at an approximate equation. The point is that the equations in continuum mechanics can be solved yet the necessary pieces may not always make it possible to achieve a solution. At any rate, this analysis alone is not enough to describe the buc kling behavior that is observed. While a buckled solution can result fr om a certain stress field being applied to a material, this tells us nothing about what causes it and the criteria necessary to make it so. PAGE 11 4 Linear Elastic Theory A polymer hydrogel is an elastic material. It has the ability to return to its initial size and shape after it has been subjected to stress. This only applies up until the point before the material yields or buckles to the applied stress. Linear elastic theory is capable of describing the behavior of a material before it yields or undergoes a lar ge deformation. If the material were to yield, its behavior then becomes nonlinear and a more complicated set of equations would be necessary to map the deformation of the material from that point on. Although this can be done, the approach typically is to initially use linear elastic theory. If this analysis fails to accurately describe the materials behavior then a nonlinear analysis is justified. The equations used in linear elastic theory are linear relationships between the stress and strain of the mat erial [57]. This along with the ability of the material to return to its initial state are the main tenets of the theory. The other tenets of the theory are that the rate at which the load is applied is of no consequence and the deformations that occur are very small provided that the material is still stable. With this added piece of knowledge, a framework can be devised that describes the behavior of a polymer hydrogel being subject to stress. The only issue is how to describe what happens to the material when it buckles without making use of the nonlinear analysis. Structural or Elastic Stability Theory Structural and Elastic Stability Theory refer to the same general idea. The idea is that there comes a point at which a material, be it elastic or not, loses its stability or its ability to support the applied load. This instability can make itself manifes t in a variety of ways. PAGE 12 5 It can be done through bending as well as through multiple buckling modes. The point at which this instability occurs depends on how much work is being done on the material as well as the amount of energy required to cause a large deformation. This means that the nature of the material and the applied forces on that material determine how it looks when it becomes unstable and what critical values of these parameters are required to make it so. This is a very crucial part in the study and analysis of buckling phenomena. This theory makes it possible to have tangible, applicable quantities that can be used to show how they affect buckling and structural stability. Biharmonic Equations and the Airy Stress Function The biharmonic equation is a differential equation that results from the equation shown below which uses the biharmonic operator 4. 4 = 0 (1 ) In Cartesian coordinates and in three dimensions the equation becomes : 4 = 2 x2 + 2 y2 + 2 z2 2 x2 + 2 y2 + 2 z2 = 4 x4 + 4 y4 + 4 z4 + 2 4 x2 y2 + 2 4 x2 z2 + 2 4 y2 z2 = 0 (2 ) There are certain complex phenomena that require the use of this equation[8]. The bending and flexing of thin flat plates is one of these phenomena. If a confined hydrogel strip is treated like a thin flat plate, then a form of the b iharmonic equation can be used to analyze the stresses acting on it. This form of the biharmonic equation is called the Airy Stress function. This function is capable of determining the equilibrium conformation of a material subjected to stress. The main idea behind this function is that at PAGE 13 6 equilibrium, the stresses acting on the object satisfy a set of compatibility criteria[5, 7] These criteria are shown below. + + + = 0 + + + = 0 + + + = 0 (3 ) With these equations in hand, the Airy Stress function can be derived. This then makes it possible for the relationships that exist within linear elastic theory to be used to further simplify the equations so that there are fewer degrees of freedom. Derivation of Airy Stress Function The Airy stress function ( ) is a functi on that describes the equilibrium of a continuum body. In the absence of body forces in two dimensions with plane strain, the strain components at equilibrium satisfy Eq. 4. 2 2 + 2 2 2 2 = 0 (4 ) From linear elastic theory the stress strain relations are: = 1 + (1 ) = 1 + (1 ) PAGE 14 7 = 1 + (5 ) Substituting the stress strain relations into Eq. 4 ( 1 ) 2 2 2 2 + ( 1 ) 2 2 2 2 2 2 = 0 (6 ) The potential functions defined by the stresses are shown below = 2 2 = 2 2 = 2 (7 ) Substituting the potential function into Eq. 6 ( 1 ) 4 4 + ( 1 ) 4 4 + 2 ( 1 ) 4 2 2 = 0 4 4 + 4 4 + 2 4 2 2 = 4 = 0 (8 ) There are other biharmonic equations that expressed in terms of strains, stresses and displacements that describe a body in equilibrium [7, 9]. It is our desire to obtain the buckled conformation of the hydrogel and this will require that the equation be expressed in terms of displacements. This same principle is applicable to the boundary conditions and they will also need to be expressed in terms of displacements so that that will be the only unknown. PAGE 15 8 Finite Difference Approximation The finite difference approximation method was used to solve the differential equation that model s the buckling of the hydrogel. This method was also used to attempt to see if changes in the parameters of the function would give rise to noticeable changes in the buckling behavior It has been implemented as a way to solve boundary value problems that c ontain partial differential equations of various types and orders [1013]. The finite difference approximation breaks down derivatives into small and gradual step changes over a given variable. For example, the first derivative of a variable y with respect to x is shown below: = ( + ) ( ) As the step size becomes smaller, the approximation of the derivative becomes more accurate. However, more calculations have to be done to achieve this accuracy. The effect of this accuracy on solution stability will be discussed in further detail in a later section. Below are the some of the finite difference approximations necessar y to solve the Airy St ress f unction. 4 4 = + 2 4 + 1 + 6 4 1 + 2 ( )4 4 4 = + 2 4 + 1+ 6 4 1+ 2( )4 4 2 2 = + 1 + 1+ + 1 1 2 + 1 2 1 ,+ 4 2 + 1 2 1+ 1 + 1+ 1 1( )2( )2 PAGE 16 9 Combining these derivative expressions make it possible to have an approximation of the Airy Stress function or any other biharmonic equation as well as the behavior of the function on the boundaries. While further derivation may be n ecessary to describe the boundary behavior, they can be accurately encompassed within the framework of finite differences. Basic Idea and S etup Once the finite difference equations have been derived, the area over which the equation applies must be defined. The points of analysis on the hydrogel are along the surface in the x and y directions. The goal is to use the overall equation to solve for the deflection of the surface in the z direction. Fig. 2 shows how the grid points are defined as well as where the boundaries lie x y Zi,jZi,j+1Zi,j1Zi+1,jZi1,j y = L y = 0 x = W x = 0 Figure 2 : Boundaries and Grid Definition of Finite Difference Method PAGE 17 10 The points are counted from the bottom left when x and y are equal to zero. They are then iteratively stepped out in increments of the defined step size upwards in the x direction. Once the width of the hydrogel is reached, the next grid point is set at when x equals zero. However, that point is one step above the first grid point in the y direction. This procedure is continued until the points mark the boundaries and the entire surface at regularly spaced intervals. Stability of Solution The solution of the finite difference approximation has multiple limitations. Since this is an approximation, there will always be a certain amount of error. There is a point however at which this error becomes intolerable to the set of equations being approximated. At this point the system can go unstable and the solution that is calculated can be very far from accurate. The von Neumann method is capable of determining the stability criteria for the finite difference equations over a specific area [10, 14]. This method introduces an error into the difference equations which is then subsequently solved for. The upper limit of this error is determined and this becomes the stability criteria. There seem to only be solutions to specific problems solved with this method and very little explanation about how to go about using the method i tself. Because of this, the number of grid points, which determine the step size of that variable, were changed and the effect on stability was analyzed. PAGE 18 11 Boundary Conditions Provided that the boundary conditions are linear in nature, a linear solving technique will provide a solution. However, non linear boundary conditions will require the use of a non linear solver. There are various ways to describe the boundaries of a material subjected to stress [6, 7, 9]. In terms of the vertical displacement in the z direction, the boundary conditions are defined as follows. For a clamped edge, = 0 (9 ) = 0 (10 ) Combining Eqs.9 & 10, = 0 (11 ) The forward finite difference approximation of the first derivative is: = + 1 ( ) (12 ) Substituting Eq. 12 into Eq.11, the expression for a clamped edge is : + 1 ( ) = 0 (13 ) For a simply supported e dge, = 0 (14 ) 2 2 + 2 2 = 0 (15 ) Combining Eqs. 14 & 15, 2 2 + 2 2 = 0 (16 ) PAGE 19 12 Expressed in forward finite differences with respect to x and in central finite differences with respect to y Eq. 16 for a simply supported edge is: + 2 2 + 1 + ( )2 + + 1 2 + 1( )2 = 0 (17 ) For a free edge, 2 2 + 2 2 = 0 (18 ) 3 3 + ( 2 ) 3 2 = 0 (19 ) Combining Eqs. 18 & 19, 3 3 + ( 2 ) 3 2 2 2 2 2 = 0 (20 ) Expressing Eq. 20 with central finite differences with respect to x and backward finite differences with respect to y the equation for a free edge is: 3 1+ 3 2+ 3( )3 + ( 2 ) + 1 2 + 1 + 1 1+ 2 1 1 1( )2( ) 2 1+ 2( )2 + 1 2 + 1 ( )2 = 0 (21 ) The boundary conditions that were used as well as where they were used are highlighted in Fig. 3. PAGE 20 13 x yy = L y = 0 x = 0 3322 322220 ZZZZ yxyyx 0 Z Z y x = W 22 220 ZZ Zxy 22 220 ZZ Z xy Figure 3 : Boundary Conditions at Specified Boundaries At the lower y boundary, the hydrogel was treated as a clamped edge while at the upper y boundary, it was treated as a free edge. Both of the x boundaries were treated as simply supported edges. It was also necessary to introduce buckling on the upper y boundary in the form of a sine curve. Without this, no solution other than z ero across the whole surface can be calculated. Zero is a very good solution to the biharmonic equation. Therefore a reason for any other solution must be justified and introduced through the relevant equations. While this induced buckling boundary condition may not accurately describe the actual buckling phenomenon, it sets up a condition which could mimic an approximation of this behavior. Matrix Inversion as a M eans of Solution Once the set of equations have been generated to map the behavior of the hydrogel strip, a solution to the conformation of the buckled surface must be obtained. Since the set of generated equations are linear in nature, a linear solving technique will be able to arrive at a solution. Matrix inversion was the method of choice due to ease of programming as well as it being a built in function in MATLAB PAGE 21 14 Matrix inversion simply requires there to be a coefficient matrix as well as a solution matrix for the defined system The coefficient matrix comes from the equations that are derived and that govern the materials behavior. The solution matrix takes into account certain known values that may exist and integrates them into finding a solution. For ex ample, if the coefficient matrix is defined as C and the solution matrix as S, the expression to solve for the variable X would be: = Hence the solution with matrix inversion is = 1 Once the solution was obtained, the values were then rearranged and plotted using the same grid point defi nition as highlighted in Fig. 2. PAGE 22 15 Analytical Solution Using the analytical solution of the biharmonic equation is a viable way to achieve numerical solutions to the buckling problem. While the results obtained through this method may not be representative of what happens experimentally, it does provide an insight into the nature of the mechanics of a very simplified case of this phenomenon. Separation of Variables The separation of variables technique was used to solve for the analytical solution of the biharmonic equation. While other methods exist, I am most familiar with this method and the solution is much easier to construct. My initial attempt at using this technique was unsuccessful After reading further, I discovered that the separation of v ariables procedure would have to be used multiple times in order to produce a solution. Otherwise, if one of the functions with respect to a certain variable is presumed, the separation is much easier. The general procedure for this functions variable separation is highlighted below. Assuming Z is a function of x and y only: ( ) = ( ) ( ) (22 ) Hence the biharmonic equation will yield PAGE 23 16 4 4 + 2 4 2 2 + 4 4 = 0 4 4 + 2 2 2 2 2 + 4 4 = 0 (23 ) Dividing by XY: 1 4 4 + 2 1 2 2 1 2 2 + 1 4 4 = 0 (24 ) New variables can be defined that can separate this function into different parts. 1( ) = 1 4 4 ; 2( ) = 1 2 2 ; 1( ) = 1 4 4 ; 2( ) = 1 2 2 (25 ) The new expression that contains these variables is: 1+ 2 22+ 1= 0 (26 ) Differentiating Eq. 26 with respect to x: 1 + 2 2 2 = 0 (27 ) Which simplifies to: 1 2 = 2 2= ( 28 ) Differentiating Eq. 26 with respect to y: 1 + 2 2 2 = 0 Which simplifies to: PAGE 24 17 1 2 = 2 2= ( 29 ) Using Eqn. 29 and defining a constant as 2, 2= 1 2 2 = 2 (30 ) Therefore: 2 2 + 2 = 0 (3 1 ) Eq. 31 has solutions of the form: ( ) = ( ) + ( ) (32 ) Using Eq. 32 and solving for Eq. 25 shows that, 1( ) = 1 4 4 = 4 (33 ) Substituting Eqs. 30 & 33 into Eq. 26 yields: 1 4 4 2 2 1 2 2 + 4= 0 (34 ) Multiply by Y: 4 4 2 22 2 + 4 = 0 (35 ) The solution to Eq. 35 can be expressed in terms of exponentials, trigonometric functions or a combination of both. Different combinations have been used and in different ways [8, 9], but the y function in terms of exponentials can be expressed as: PAGE 25 18 ( ) = 1+ 2 + 3 + 4 (36 ) Therefore, the overall expression for Z u sing Eqs. 32 & 36 is : ( ) = [ ( ) + ( ) ] 1+ 2 + 3 + 4 (37 ) From observed experiments it is fair to assume that the buckling modes take on a sinusoidal shape as a function of x. Hence the known function behav ior will be a sine function with respect to x. Implementing this into the biharmonic equation in terms of displacement yields the following results. 4Z = 4Z x4 + 4Z y4 + 2 4Z x2 y2 = 0 Z = ( ) ( ) ; ( ) = sin ( ) 4Z = 4sin ( ) ( ) + sin ( ) 4 ( ) y4 2 2sin ( ) 2 ( ) y2 = sin ( ) 4 ( ) y4 2 22 ( ) y2 + 4 ( ) = 0 (38 ) The final result shows that the variables can be separate. However, we still do not know w hat the y function is. If Eq. 38 is divided by sin( mx), the result is: 4 ( ) y4 2 22 ( ) y2 + 4 ( ) = 0 (39 ) Solving Eq. 39 will give rise to a similar solution as Eq. 35. Using one of the possible solutions, the y function can be expressed as: ( ) = 1+ 2 + 3 + 4 (40 ) The overall expres sion for Z would then be: = sin ( ) [ 1+ 2 + 3 + 4 ] (41 ) PAGE 26 19 Now that the there is an expression for Z, exact values have to be determined for the constants. This procedure is detailed in the following section. Using Boundary C onditions to S olve for Constants There are known solutions to the biharmonic equation which can be used directly. The previous two sections highlighted methods as to how these solutions can be obtained. The next step requires the use of boundary conditions to solve for the constants followed by the functions implementation. The series solution, shown below, is another example of an existing solution that can be used to solve this problem. ( ) = = 0 = 0 (42) Where Then only thing that would be needed to use this solution are values of Cmn. For the case of a uniaxial compression in the x direction, a modified form of the biharmonic equation is required. This function is shown below. 4 = 1 2 2 + 2 2 + 2 2 (43 ) Since the only force acting is in the x direction, NY and NXY are set equal to zero. Also, because the force is a compression force, NX should undergo a sign change. Hence the expression becomes: 4 = 2 2 (44 ) Due to this modification, the y function will take on a different form. PAGE 27 20 ( ) = 1 + 2+ 3 ( ) + 4 ( ) (45 ) = ; = 2+ 2 ; = 2+ 2 Hence the expression for Z will be: = sin ( ) [ 1 + 2+ 3 + 4 ] (46 ) The boundary conditions from Eqs. 9, 10, 18 &19 provide four independent equations for the four unknown constants. These boundary conditions as well as the expressions they provide are shown below. When = 0 = 0 1+ 2+ 3= 0 (47 ) When = 0 = 0 1+ 2+ 4= 0 (48 ) When = ,2 2 + 2 2 = 0 ( 2 2 ) 1+ ( 2 2 ) 2+ ( 2cos ( ) 2 cos ( ) ) 3+ ( 2sin ( ) 2 sin ( )) 4= 0 (49 ) When = ,3 3 + ( 2 )3 2 = 0 ( 3 + ( 2 ) 2 ) 1+ 3 ( 2 ) 2 2+ 3sin ( ) + ( 2 ) 2 sin ( ) 3+ 3cos ( ) ( 2 ) 2 cos ( ) 4= 0 (50 ) PAGE 28 21 These equations are linear and hence can be solved using a linear solver. The coefficient matrix used for this set of equations is shown below. 22 22 3232 22 22 323211 22 10 0 coscos sinsin sin2sincos2cosLL LL mm LLLL mm mm mmee ee eeee LLLL LLLL T he solution to a buckling problem requires that the determinant of the coefficient matrix for this set of equations to have a determinant of zero[2, 4, 9]. This means that the matrix cannot be inverted since it will be a singular matrix. Therefore matrix inversion will not be capable of arriving at a solution. However, the matrix can be modified to reduced row echelon form which will then leave four simplified equations. The fourth equation does not pr ovide any relevant information so there are only three useful equations remaining. Because of this, any nonzero value for the fourth constant can be used and the rest of the equations can then be satisfied. The other three constants will be factors of the fourth constant because of this. The expression for the determinant of this system is: 2232 3222 223 2 2232coscoscos2cos sin2sinsinsin sin2sin coscos 2mm mm LL mm LL mmLLLL LLLL eeLL LLee PAGE 29 22 2232 3222 223 2 2232coscoscos2cos sin2sinsinsin sin2sin coscos 2mm mm LL mm LL mmLLLL LLLL eeLL LLee 223 2 2232 223 2 2232 2232 2cos2cos sinsin 2 cos2cos sinsin 2 2LL mm LL mm LL mm LL mm LLL L mmeeLL LLee eeLL LLee eeee e 232(51) 2LLL L mmeee Once values have been selected that give rise to a determ inant equal to zero, using Eq. 51, the value of the constants in the buckling equation can be obtained. With the function well defined, an analysis can now be done on the variables that are present within the equation. This analysis will give rise to critical values and show where the threshold of stability lies due to a change in certain function parameter values. PAGE 30 23 Determination of Critical Buckling Values A key aspect of buckling analysis is the amount of stress a body can endure before is goes unstable and buckles. This analysis provides an insight into the effect that certain material properties have on the material s stability. By varying properties like the hydrogels length, width, modulus of rigidity, P oisson ratio, and magnitude of compression, the effect of these changes on buckling can be observed. Likewise, the critical values of these parameters can also be determined. The critical value of a parameter represents a point where an increase or decrease from this point will cause instability. Whether that change should be an increase or a decrease depends on whether the plot of the parameter is a minimum or maximum curve. The methods required to obtain these values are detailed in the following sections. Failure of Plane Strain Assumption The plane strain assumption implies that there is no reason for there to be displacement in any direction other than those within the plane being subjected to stress [5, 7]. For example, if a material is being stressed in the x and y directions, then there should only be displacement in those directions. However, if the magnitude of the applied force increases, there is a point at which the plane strain assumption ceases to apply. This is when buckling occurs This is due to an instability in the solution of the biharmonic equation which occurs when the determinant of the coefficient matrix for the equations constants equals zero. From this particular moment, a lot of information can be PAGE 31 24 determined. Various parameters within the equation can be changed and their effect on the buckling solution can be observed. Determinant Method for Critical Buckling Values Using the analytical solution to the biharmonic equation, the magnitude of the critical force or pressure required for buckling can be determined. Similar work has been done which verifies that these quantities can be obtained[2]. Eq. 51 is the expression that represents the determinant for the solution to the constants in Eq. 46. As mentioned before, setting this determinant to zero will be representative of the instability and the buckling phenomenon. Once this is done, it is desired to determine how much compression force is necessary to cause different buckling modes to occur. So, using the same determinant expression, the value of m was changed. Simultaneously, with each change of m, the minimum value of NX was solved for. It is crucial that the value of NX be the lowest possible sol ution to Eq. 51. Fig. 4 shows that there are multiple solutions that exist that will result in a determinant equal to zero. Since it is the goal to obtain the critical value of NX, we have to observe what the smallest possible value is that will cause buckling to occur. PAGE 32 25 Figure 4 : Determinant Behavior w.r.t. Nx The result of this procedure gives rise to a pl ot like the one shown in Fig. 5. Clearly there is a minimum value of NX and a corresponding value of m to go with it. This minimum point signifies the minimum amount of compressive force required to buckle. While the value of m at that point represents the number of half waves that would be expected in the buckled conformation. PAGE 33 26 Figure 5 : Plot of Solution of Critical Values of Nx This procedur e can be further expanded to determine the critical values of other parameters. However, just like the compression force, the determinant must be analyzed so that the minimum value of the parameter is used and not just any value that satisfies the determinant criteria. The finite difference method could also be used to determine the critical values. A nonzero solution to the deflection surface would also come about when the determinant of the coefficient matrix is equal to zero. The only difference is the number of equations that will have to be used in order to solve for the values. For the sake of simplicity and ease of programming, this procedure was not used. PAGE 34 27 Energy Method The forces of compression acting on the hydrogel perform work on the hydrogel. Beca use of this there is a certain amount of work that once exceeded, will give rise to an instability and cause the onset of buckling. This energy is called the strain energy It is defined below with the following equation[9]: = 2 2 2 + 2 2 2 2 (1 ) 2 2 2 2 2 2 (52 ) There will then be a minimum value of the strain energy that will give rise to a critical amount of buckling. W hat this amount is depends on how much work is done. The equation for the work done on the hydrogel is shown below: = 1 2 2+ 2+ 2 (53 ) Since there is only compression in the x direction = 1 2 2 (54 ) To determine what the critical magnitude of compression is, both of the work and the strain energy have to be equal at a minimum value of the strain energy. So setting the work and strain energy equal to one another: = PAGE 35 28 2 2 2 + 2 2 2 2 (1 ) 2 2 2 2 2 2 = 1 2 2 (55 ) Therefore: 2 2 2 + 2 2 2 2 ( 1 ) 2 2 2 2 2 2 + 1 2 2 = 0 (56 ) The solutions to the work and strain energy equations are very complicated. The solutions were obtained using the Symbolic Math Toolbox in M ATLAB. Since an analytical solution exi sts for Z, this was used in the equations and the integrals as well as the derivatives were then evaluated. The integrals were also evaluated within their limits of zero to W and zero to L for x and y respectively. Eq. 56 can be used to solve for the necessary critical values in the same manner as with Eq. 51. The value of m can be changed and the desired critical value can be solved for. Fig. 6 shows a plot of the work and strain energy coming together at a particular point. That point is the critical val ue of NX when m is one and this value is in agreement with the value calculated through the determinant method. PAGE 36 29 Figure 6 : Work and Strain Energy Plots Used for Critical Values of Nx This method can also be used to determine the critical values of other parameters. D ue to the length and complexity of the programming code, the determinant method was the preferred method of solution. PAGE 37 30 Bending and Wasted Modes of Buckling When a swollen hydrogel buckles it has been observed that this can occur in one of two ways [2, 15]. These types of buckling can be called bending and wasted modes of buckling. Both of these modes are shown in Fig. 7 Figure 7 : a ) Wasted and b) Bending Modes of Buckling The reason for this difference in buckling behavior is not well understood. Hence it is desired to possibly find out w hy this difference occurs as well as what factors affect it. Equations I n volved with Different Buckling Modes An analysis exists which can describe the bending and wasted buckling modes. This analysis is similar to what was done with the analytical solution, however there are some differences. This analysis defines the flat plate or strip being in the x and z pl ane. The compression is still applied in the x direction and all other characteristics are the same. PAGE 38 31 Looking at the problem from this perspective will show a view from the top of how the hydrogel buckles. The equations required for this analysis are shown below The equations used in the biharmonic equation for a compressible and incompressible strip respectively are: = ( x ) ( ) ( 57 ) = ( x ) ( ) (58 ) The overall governing equation that applies to compressible and incompressible strips for the z function is: 2 2 ( 2)2 2 2 ( 3)2 ( ) = 0 4 ( ) 4 ( 2)22 ( ) 2 ( 3)22 ( ) 2 + ( 2)2( 3)2 ( ) = 0 ( 59 ) If 23 then the equations for the bending and wasted modes of buckling respectively are: ( ) = 1 ( 2) + 2 ( 3) (60 ) ( ) = 1 ( 2) + 2 ( 3) (61 ) If 23 then the equations for the bending and wasted modes of buckling respectively are: ( ) = 1 ( 2) + 2 ( 3) (62 ) ( ) = 1 ( 2) + 2 ( 3) (63 ) PAGE 39 32 Similar to what was highlighted before, the method for solving the constants in this equation require the use of the boundary conditions. With the boundary conditions expressed, a set of linear equations exist which can hence be solved. The instability results when the determinant of the coefficient matrix for this problem equal zero. Unfortunately this analysis tells us nothing about why these phenomena occur. All that is presented are equations that mimic both phenomena separately. In my opinion, one equation should be capable of showing how there could be a transition from one type of buckling to the other based on certain property changes within the hydrogel strip. Possible Reasons for Different Buckling Modes It has also been observed that the dimensions of the confined hydrogel play a distinct role in determini ng which of these modes occur [2] When the hydrogels thickness in the z direction increases, there is an overall decrease in the amount of buckling that is observed. This could be due to a number of reasons but one of them seems to be very reasonable. Figure 8 : Buckling Transition from One Mode to Another Observe the hydrogel in Fig. 8. As the hydrogel is thickened, its volume is increased assuming that it is still the same height. Yet since it is still a part of the same hydrogel, PAGE 40 33 when it is swollen, it undergoes the same magnitude of compression as the rest of the gel. However, the part of the gel with a greater thickness will undergo les s pressure while the thinner part will undergo more. The greatest magnitude of buckling is observed at the free edge which is farthest away from the confined surface. This causes a gradient in buckling which would then be a function of the thickness or the amount of pressure that that region experiences uniaxially. T he magnitude of compression in the z direction may also increase as the hydrogel gets thicker. To analyze this theory a more complicated analysis is necessary which will require three dimensional analysis of the buckling problem. While the theoretical framework exists for three dimensional analysis of buckling, substantial modifications would have to be made to fit the theory to this problem. PAGE 41 34 Results With all of the mentioned methods of solving the buckling problem, the goal is to determine approximately how much buckling can be expected for a specific size of confined polymer. Also, noting how changing different properties of the hydr ogel can give ris e to different magnitudes or types of buckling. Finite Difference Solution The finite difference solution is only capable of arriving at an approximation of the buckled defl ection surface of the hydrogel. After using the appropriate boundary conditions t o describe the gel, the approximation seems to only be able to mimic the edge undulation or the wasted form of the buckling. Curvature was induced on the free edge and the boundary conditions described in Eqs. 9,10,14,15,18 and 19 were used as well. Combinations of forward, backward, and central finite difference schemes were used to define the coefficient matrix. Trial and error methods were used to determine the combinations of finite differences that made the system most stable. A spy of the coefficient matrix is shown in Fig. 9. The dots represent the places in the matrix where there are values. This coefficient matrix is solving for a total of one hundred grid points arranged in a five by twenty grid. PAGE 42 35 Figure 9 : Spy of the Coefficient Matrix of Finite Difference Problem Due to the fact that the stability of the finite difference solution is of great importance, limitations exist in regards to the number of grid points that can be used. Fig. 10 shows how changes in the number of grid points in the y direction affect the solution. The grid points were increased from four, to five, and then to six. When there are more than five grid points, the solution is no longer stable. PAGE 43 36 Figure 10: Solution Instabilities with Different Numbers of Grid Points on the y Axis. Twenty grid points were used on the x Axis and a) four b) five and c) six points were used on the y axis. Fig. 11 shows how the number of grid points on the x direction affect the stability of the solution. The grid points were increased from twenty, to thirty, and then to thirty four. Thirty grid points in this direction seem appropriate. While thirty two grid points seem to be the threshold of stability, there is no need to perform the unnecessary calculations as no new information is gained by doing this. PAGE 44 37 Figure 11: Solution Instabilities with Different Numbers of Grid Points on the x Axis. Five grid points were used on the y Axis and a) twenty b) thirty and c) thirty four were used on the x axis. Fig. 11b shows the most stable solution to the problem with the greatest amount of grid points allowed. There are thirty grid points in the x direction and five in the y direction. Analytical Solution The analytical solution was used to perform multiple tasks. It was first used to verify the buckled deflection surface predicted by the finite difference approximation. Then, the critical values of the wavelength and compressive stress were determined so that relationships between the parameters in the equation could be found. Fig. 12 shows the analytical solution of the buckled deflection surface. PAGE 45 38 Figure 12: An Analytical Solution with Five Half Waves Upon changing the parameters in the analytical equation, there was no noticeable change in the type of buckling that was observed. For this reason, the analytical solution as it is seems to only be capable of mimicking the bending mode of buckling just like the finite difference solution Determinants The analysis of the behavior of the determinants with respect to each parameter is very necessary. If changes in the parameter do not cause the determinant to equal zero, then that parameter cannot cause buckling to occur. This can also show ahead of time what can be expected for the critical buckling values. Fig. 13 shows a plot of how the determinant behaves using the determinant method to solve for the critical values of Nx. There are two nonzero solutions that satisfy the determinant criteria for every value of m. The minimum value of the solution was selected since this value would represent the minimum possible amount of compression required for buckling. PAGE 46 39 Figure 13: Determinant Behavior w.r.t. Nx at Different Values of m Fig. 14 shows the effect that the flexural rigidity has on the determinant. Unlike any of the other parameter analyses, the flexural rigidity has only one solution that satisfies the determinant criteria. PAGE 47 40 Figure 14: Determinant Behavior w.r.t. Flexural Rigidity at Different Values of m Fig. 15 shows how the ratio of length to width affects the determinant. Upon further analysis, it was observed that there are infinite solutions that satisfy the determinant criteria with the length and w idth ratio. As the length to width ratio is increased, the behavior of the determinant becomes more and more unstable. Increasing the value of the length to width ratio causes the values of the determinant to deviate further to higher amplitudes in an unstable sine curve way. PAGE 48 41 Figure 15: Determinant Behavior w.r.t. Length and Width Ratio Analysis of the Poisson ratio shows that this parameter can have an effect on the buckling behavior. However the points at which the solutions occur are out of range for an elastic material. An elastic materials Poisson ratio i s between zero and one and there were very few solutions that fell within this range. Critical Values All of the determinant values for the respective parameters highlighted in the previous section were solved for at different values of the wavenumber. The critical value plots all have a unique characteristic. There is either a minimum or a maximum point w hich highlights a limit for that parameter. What that limit means or represents depends on the parameter being analyzed. PAGE 49 42 Fig. 16 shows the critical values of flexural rigidity obtained at different values of Nx. Figure 16: Critical Values of Flexural Rigidity at Different Values of m and Nx All of the critical flexural rigidity curves have a maximum value. This maximum is the largest possible value of flexural rigidity that can cause buckling. This, in theory, does make sense. If an object is more rigid, it is less likely to buckle. So there should be an upper limit of this value, above which no buckling will occur. Increasing the magnitude of compression increases this maximum value. This also makes sense because a larger f orce of compression will be able to counter the effects of a more rigid body or a greater value of flexural rigidity. Fig. 17 shows the critical values of Nx obtained at different values of flexural rigidity. PAGE 50 43 Figure 17: Critical Values of Nx at Different Values of m and Flexural Rigidity Si milar to the analysis of Fig. 16, the behavior of the critical values can be explained. An increase in the value of flexural rigidity will result in a greater force of compression necessary to cause buckling. In both Figs. 16 & 17, these changes have no effect on the wavenumber at which buckling occurs. Hence these changes have no effect on the wavelength or the number of observed buckles on the deflection surface. Fig. 18 shows the critical values of Nx and how they are affected by changes in the length to width ratio. PAGE 51 44 Figure 18: Critical Values of Nx at Different Values of m and Length to Width Ratio Assuming that the length of the hydrogel is kept at a constant value, the width of the hydrogel can be changed to give rise to different values of the length to w idth ratio. So an increase in the width of the hydrogel will give rise to a decrease in the length to width ratio. Therefore, Fig. 18 shows that increasing the width of the hydrogel causes an increase in the minimum critical amount of compression necessary for buckling. The changes in the length to width ratio also cause a shift in the critical buckling wavenumber. Increasing the width of the hydrogel causes the critical wavenumber to increase. This means that increasing the width of the hydrogel will cause an increase in the number of observed buckles on the deflection surface. Fig. 19 shows the critical values of the length to width ratio and how changes in Nx affect them. PAGE 52 45 Figure 19: Critical Values of Length to Width Ratio at Different Values of m and Nx This result is in direct correspondence to res ults obtained in Fig. 18. A wider hydrogel requires a greater amount of compression to cause buckling. Therefore more compression will cause a decrease in the critical ratio of length to width. The increase in compression will also cause an increase in the critical buckling w avenumber This is equivalent to the relationship in Fig. 18 where a decrease in the length to width ratio also causes an inc rease in the critical buckling wavenumber. The obtained results overall, highlight key relationships that exist between the parameters that are present in the thin plate buckling analysis. In reality however, these parameters may not necessarily be so separate and easy to control. If they can be controlled, the general trends of the observed relationships should still retain their validity. PAGE 53 46 Conclusion The buckling phenomenon, in reality, is quite complex. However, assumptions can be made to simplify the problem and make analysis more feasible. Treating a surface confined polymer hydrogel like a thin plate subject to compression is a way to simplify the problem yet find out something about its buckling behavior. There are different methods of solving the equations governing this problem, but the desire to obtain quantities that are representative of this phenomenon requires the use of the analytical solution to the problem. The relationships obtained from the analytical solution were able to show how the flexural rigidity, magnitude of compression, and the length to width ratio affected the buckling behavior. The length to width ratio is a parameter that i s easy to control. However, being able to change the flexural rigidity and the magnitude of compression may introduce a new set of factors that would have to be analyzed. For example, increasing the cross link density within the hydrogel may cause an incre ase in the value of the flexural rigidity. But this increase in cross link density may also limit the amount of water that can diffuse into the hydrogel. The amount of diffused water could very well determine the amount of compression that the hydrogel undergoes. As a result, both of these parameters, flexural rigidity and the magnitude of compression, depend on the cross link density. The relationship that exists between these quantities will have to be determined in order to have the desired control over them and hence the desired amount of buckling. Other relationships that are coupled to PAGE 54 47 a single parameter may exist and they must also be analyzed to see how they affect the key buckling parameters. The current framework does not however solve for the ampl itude of the buckled solution. This is an important factor that is necessary to describe buckling behavior. But this parameter is definitely a function of the swelling ratio of the gel before and after swelling. If the gel swells to a size substantially gr eater than its initial size, then the amplitude of the buckles should increase. Experiments would have to be done to determine what the nature of this relationship may be as it is not included in the thin plate assumption. Another key analysis is the descr iption of the bending and wasted modes of buckling. As was mentioned before, equations exist that describe the two phenomena separately. But a method of analysis must exist that shows how the transition would occur from one mode to the other. From what I have observed, the magnitude of compression in the z direction increases as the hydrogel gets thicker and this could explain the two different types of buckling. However, the appropriate equations would have to be used in order to verify that this idea has validity. Due to a lack of time, further analysis in this regard could not be pursued. Further experiments may also have to be done to verify that the described relationships from the thin plate assumption are accurate. For example, an experiment where the length of the gel is kept constant, yet different widths of the hydrogel are used would help to verify the length to width ratio relationships. PAGE 55 48 References 1. Galaev, I.Y. and B. Mattiasson, 'Smart' polymers and what they could do in biotechnology and medicine. Trends in Biotechnology, 1999. 17(8): p. 335340. 2. Mora, T. and A. Boudaoud, Buckling of swelling gels. European Physical Journal E, 2006. 20(2): p. 119124. 3. Russell, T.P., Surfaceresponsive materials. Science, 2002. 297(5583): p. 964967. 4. Chajes, A., Principles of structural stability theory Civil engineering and engineering mechanics series. 1974, Englewood Cliffs, N.J.,: PrenticeHall. xii, 336 p. 5. Lai, W.M., D. Rubin, and E. Krempl, Introduction to continuum mechanics 3rd ed. 1993, Oxford ; New York: Pergamon Press. xiv, 556 p. 6. Landau, L.D. and E.M. Lifshitz Theory of elasticity Their Course of theoretical physics,. 1959, London, Reading, Mass.: Pergamon Press; AddisonWesley Pub. Co. 134 p. 7. Saada, A.S., Elasticity theory and applications 2nd ed. 1993, Malabar, Fla.: Krieger. xvi, 775 p. 8. Selvadurai, A.P.S., Partial differential equations in mechanics 2000, Berlin ; New York: Springer. 9. Timoshenko, S., Theory of elastic stability 2d ed. Engineering societi es monographs. 1961, New York,: McGraw Hill. 541 p. 10. Constantinides, A. and N. Mostoufi, Numerical methods for chemical engineers with MATLAB applications Prentice Hall international series in the physical and chemical engineering sciences. 1999, Upper Saddle River, N.J.: Prentice Hall PTR. xxv, 559 p. 11. Altas, I., J. Erhel, and N.M. Gupta, High accuracy solution of threedimensional biharmonic equations. Numerical Algorithms, 2002. 29(13): p. 119. 12. Buzbee, B.L. and F.W. Dorr, Direct Solution of Biharmonic Equation on Rectangular Regions and Poisson Equation on Irregular Regions. Siam Journal on Numerical Analysis, 1974. 11(4): p. 753763. PAGE 56 49 13. Ehrlich, L.W. and M.M. Gupta, Some Difference Schemes for Biharmonic Equation. Siam Journal on Numerical Analysis, 1975. 12(5): p. 773790. 14. Kelley, W.G. and A.C. Peterson, Difference equations : an introduction with applications 2nd ed. 2001, San Diego: Harcourt/Academic Press. ix, 403 p. 15. Guz*, A.N., Fundamentals of the threedimensional theory of stability of deformable bodies Foundation of engineering mechanics. 1999, Berlin ; New York: Springer. xvi, 555 p. PAGE 57 50 Appendices PAGE 58 51 Appendix A: Nomenclature An, Cn Cmn Equation Constant C Coefficient Matrix D Flexural Rigidity exx ,yy,xy Strain Tensor Component Fx,y,z Body Force L Length of Hydrogel in y D irection m Wavenumber (Number of Half Nx Uniax ial Compression Force in the x D irection S Solution Matrix W Width of hydrogel in x D irection Z Solution to Buckled Deflection Surface m, Equation Parameters (defined in the next section) Poisson Ratio Equation Parameter xx,yy,xy Stress Tensor Component Potential Function Equation Constant PAGE 59 52 Appendix B: Derivations of Finite Difference Equations Starting off with the finite difference equation of a second derivative: 2 2 = + 1 2 + 1 ( )2 The third order derivative is determined via the following equations : 3 3 = 2 2 + 2 2 2 2 3 3 = 2 2 + 2 2 2 2 Expressions exist for the above third order derivatives. The following derivative was derived as shown below. 3 2 = 2 2 + 2 2 2 2 = + 1 + 1 2 2 + 1 2 + 1 + 1 2 + 1 1 2 2 1 2 + 1 1 2 ( )2( ) The fourth order derivative is determined via the following equations : 4 4 = 3 3 + 2 3 3 2 4 4 = 3 3 + 2 3 3 2 Expressions already exist for the first two derivatives. However the following derivative was derived as shown below. PAGE 60 53 Appendix B Continued 4 2 2 = 3 2 + 2 3 2 2 = + 1 + 1 2 + 1+ 1 + 1 + 1 + 2 1 ( )2( )2 + 1 2 + 1 + 1 1+ 2 1 1 1 ( )2( )2 = + 1 + 1 2 + 1+ 1 + 1 2 + 1 + 4 2 1 + + 1 1 2 1+ 1 1( )2( )2 PAGE 61 54 Appendix C: Derivation of Biharmonic Equation in Terms of Displacement The shearing forces are defined as: = = The projection of the forces in the z direction are: 2 2 + 2 2 + 2 2 + + The overall equilibrium expression including the addition of all the projected forces and a load (q) is: + + + 2 2 + + 2 2 + + 2 2 + + = 0 Or, + + + 2 2 + + 2 2 + + 2 2 + + = 0 Since Nx, Ny, and Nxy are constants, the expression becomes: + = + 2 2 + 2 2 + 2 2 Us ing the expression for the shearing forces, 2 2 2 + 2 2 2 = + 2 2 + 2 2 + 2 2 PAGE 62 55 Appendix C Continued 2 2 2 2 + 2 2 = + 2 2 + 2 2 + 2 2 The moments are defined as: = 2 2 + 2 2 = 2 2 + 2 2 = (1 ) 2 Substituting the moments into the expression: 2 2 2 2 + 2 2 + 2 (1 ) 2 2 + 2 2 2 2 + 2 2 = + 2 2 + 2 2 + 2 2 4 4 + 4 2 2 + 2 4 2 2 2 4 2 2 + 4 4 + 4 2 2 = 1 + 2 2 + 2 2 + 2 2 4 4 + 2 4 2 2 + 4 4 = 1 + 2 2 + 2 2 + 2 2 In the absence of a load and compression in directions other than the x direction, the expression simplifies to: 4 = 1 2 2 If the force acting is a compression, the sign of Nx will be negative. Otherwise there is no need for a sign change. PAGE 63 56 Appendix D: Proof of Solution to Biharmonic Equation If the proposed solution is as follows = sin ( ) [ 1 + 2+ 3 + 4 ] = ; = 2+ 2 ; = 2+ 2 The biharmonic equation with a compression force in the x direction is: 4Z = 4Z x4 + 4Z y4 + 2 4Z x2 y2 = NxD 2Z x2 Solving for the LHS and dividing by sin( mx): 4[ 1 + 2+ 3 + 4 ] 2 2[ 21 + 22 23 24 ] + [ 41 + 42+ 43 + 44 ] Which further simplifies to: ( 4 2 22+ 4) 1 + ( 4 2 22+ 4) 2 + ( 4+ 2 22+ 4) 3 + ( 4+ 2 22+ 4) 4 Evaluating the coefficients: ( 4 2 22+ 4) = 4 2 2 2+ 2 + 2+ 2 2= 4 2 4 2 2 2 + 4+ 2 2 2 + 2= 2 PAGE 64 57 Appendix D Continued ( 4+ 2 22+ 4) = 4+ 2 2 2+ 2 + 2+ 2 2= 4 2 4+ 2 2 2 + 4 2 2 2 + 2= 2 Hence the LHS is: 2[ 1 + 2+ 3 + 4 ] The RHS is 2 2 = [ 2sin ( ) ][ 1 + 2+ 3 + 4 ] Dividing by sin( mx) just as with the LHS, RHS then becomes: 2[ 1 + 2+ 3 + 4 ] Since the left and right hand sides are equal, the proposed solution satisfies the bi harmonic equation. PAGE 65 58 Appendix E: Derivation of the Strain Energy Equation The expression for the strain energy is obtained by analyzing the contributions due to both bending and twisting. The expression for the strain energy due to bending is: = 1 2 2 2 + 2 2 The bending moments with respect to x and y are defined as: = 2 2 + 2 2 = 2 2 + 2 2 Substituting the moments into the strain energy equation due to bending, = 1 2 2 2 + 2 2 2 2 2 2 + 2 2 2 2 = 2 2 2 2+ 2 2 2 2 + 2 2 2+ 2 2 2 2 = 2 2 2 2+ 2 2 2 2 2 + 2 2 2 The strain energy due to twisting is defined as: = 2 The twisting moment is defined as: = (1 ) 2 Substituting the twisting moment into the strain energy due to twisting equation, = (1 ) 2 2 Combining the strain energies due to bending and twisting, PAGE 66 59 Appendix E Continued = + = 2 2 2 2+ 2 2 2 2 2 + 2 2 2 + (1 ) 2 2 = 2 2 2 2+ 2 2 2 2 2 + 2 2 2+ 2 (1 ) 2 2 The overall strain energy is obtained by integrating the overall strain energy expression. = 2 2 2 2+ 2 2 2 2 2 + 2 2 2+ 2 (1 ) 2 2 According to Timoshenko, the overall expression can be expressed as [reference]: = 2 2 2 + 2 2 2 2 ( 1 ) 2 2 2 2 2 2 This expression is not algebraically equivalent to the other expression yet it still yields solutions that are in agreement with the determinant method. PAGE 67 60 Appendix F: Determinant of a 4x4 Matrix If a 4x4 matrix is defined as: 11121314 21222324 31323334 41424344aaaa aaaaaaaa aaaa The determinant will be 222324 212324 212224 113233341231333413313234 424344 414344 414244 212223 14313233 414243det det det det aaaaaaaaa aaaaaaaaaaaa aaaaaaaaa aaa aaaa aaa The determinants of the 3x3 matrices are as follows: 222324 323334223344344323324434422432433342 424344det aaa aaaaaaaaaaaaaaaaaa aaa 212324 313334213344344323314434412431433341 414344det aaa aaaaaaaaaaaaaaaaaa aaa 212224 313234213244344222314434412431423241 414244det aaa aaaaaaaaaaaaaaaaaa aaa PAGE 68 61 Appendix F Continued 212223 313233213243334222314333412331423241 414243det aaa aaaaaaaaaaaaaaaaaa aaa Hence the overall expression for the determinant is: 11223344344323324434422432433342 12213344344323314434412431433341 13213244344222314434412431423241 1421324333422231aaaaaaaaaaaaaaaa aaaaaaaaaaaaaaaa aaaaaaaaaaaaaaaa aaaaaaaa 4333412331423241aaaaaaaa PAGE 69 62 Appendix G: MATLAB Code Finite Difference Approximation of Biharmonic Equation function Buckling_2D % CENTRAL FINITE DIFFERENCE APPROXIMATION % OF THE BIHARMONIC EQUATION % % by: Abiola Shitta % (Grid defined below) % % Z(i,j+2) % Z(i1,j+1) Z(i,j+1) Z(i+1,j+1) %Z(i2,j) Z(i1,j) Z(i,j) Z(i+1,j) Z(i+2,j) % Z(i1,j1) Z(i,j1) Z(i+1,j1) % Z(i,j2) %Used to solve: % % % d4Z d4Z d4Z Nx d2Z % + + 2+ = 0 % dx4 dy4 dx2dy2 D dx2 clc; clear; close all; %Width of Gel (x direction)% W=30; %Length of Gel (y direction)% L=5; %Thickness of Gel (z direction)% H=0.5; %No of Grid Points along Width% NW=20; %No of Grid Points along Length% NL=5; %Step Size along Width% dw=W/(NW1); %Step Size along Length% dl=L/(NL1); Nx=100; %Poisson Ratio pr=0.25; %Modulus of Elasticity E=10*(12*(1pr^2))/(H^3); %Flexural Rigidity D=E*H^3/(12*(1pr^2)); N=NL*NW; PAGE 70 63 Appendix G Continued Z=zeros(N); %Coefficients% a=dw^4; b=2*dw^2*dl^2; c=4*dw^44*dw^2*dl^2; d=dl^4; e=4*dl^44*dw^2*dl^2+(Nx/D)*dw^2*dl^4; f=6*dl^4+6*dw^4+8*dw^2*dl^42*(Nx/D)*dw^2*dl^4; %Initial Coefficient Matrix for i=1:N if i+2*NW<=N Z(i,i+2*(NW))=a; %Z(i,j+2) end if i+1+NW<=N Z(i,i+1+(NW))=b; %Z(i+1,j+1) end if i+NW<=N Z(i,i+1*(NW))=c; %Z(i,j+1) end if i+1NW>0 Z(i,i+1(NW))=b; %Z(i+1,j1) end if i 2>0 Z(i,i2)=d; %Z(i2,j) end if i 1>0 Z(i,i1)=e; %Z(i1,j) end Z(i,i)=f; %Z(i,j) if i+1<=N Z(i,i+1)=e; %Z(i+1,j) end if i+2<=N Z(i,i+2)=d; %Z(i+2,j) end if i 1+(NW)<=N Z(i,i1+(NW))=b; %Z(i1,j+1) end if i 1*NW>0 Z(i,i1*(NW))=c; %Z(i,j1) end if i 1 NW>0 Z(i,i1 (NW))=b; %Z(i1,j1) end if i 2*NW>0 Z(i,i2*(NW))=a; %Z(i,j2) end end %Initial Solution Matrix% C=zeros(N,1); %BOUNDARY CONDITIONS% %Lower y boundary (dZ/dyZ=0) (FFD) for i=2:NW1 PAGE 71 64 Appendix G Continued Z(i,:)=0; if i+NW<=N Z(i,i+1*(NW))=1; %Z(i,j+1) end Z(i,i)=1 dl; %Z(i,j) end %Lower x boundary (FFD) for i=1:NL Z(NW*(i1)+1,:)=0; if NW*(i1)+1+2<=N Z(NW*(i1)+1,NW*(i1)+1+2)=dl^2; %Z(i+2,j) end if NW*(i1)+1+NW<=N Z(NW*(i1)+1,NW*(i1)+1+1*(NW))=pr*dw^2; %Z(i,j+1) end Z(NW*(i1)+1,NW*(i1)+1)=dl^22*pr*dw^2 dw^2*dl^2; %Z(i,j) if NW*(i1)+1+1<=N Z(NW*(i1)+1,NW*(i1)+1+1)=2*dl^2; %Z(i+1,j) end if NW*(i1)+11*NW>0 Z(NW*(i1)+1,NW*(i1)+11*(NW))=pr*dw^2; %Z(i,j1) end end %Upper x boundary (BFD) for i=1:NL Z(NW*i,:)=0; if NW*i+NW<=N Z(NW*i,NW*i+1*(NW))=pr*dw^2; %Z(i,j+1) end if NW*i1>0 Z(NW*i,NW*i1)=2*dl^2; %Z(i1,j) end if NW*i2>0 Z(NW*i,NW*i2)=dl^2; %Z(i2 ,j) end Z(NW*i,NW*i)=dl^22*pr*dw^2 dw^2*dl^2; %Z(i,j) if NW*i1*NW>0 Z(NW*i,NW*i1*(NW))=pr*dw^2; %Z(i,j1) end end %Upper y boundary (CFD) for i=(NL1)*NW+2:N1 Z(i,:)=0; if i+2*NW<=N Z(i,i+2*(NW))=dw^2; %Z(i,j+2) end if i+NW<=N Z(i,i+1*(NW))=2*dw^22*(2pr)*dl^22*dw^2*dl; %Z(i,j+1) end if i 1 NW>0 Z(i,i1 (NW))=(2pr)*dl^2; %Z(i1,j1) end if i 1+(NW)<=N Z(i,i1+(NW))=(2pr)*dl^2; %Z(i1,j+1) end PAGE 72 65 Appendix G Continued if i 1>0 Z(i,i1)=2*pr*dl^3; %Z(i1,j) end Z(i,i)=4*dw^2*dl4*pr*dl^3; %Z(i,j) if i+1<=N Z(i,i+1)=2*pr*dl^3; %Z(i+1,j) end if i+1+NW<=N Z(i,i+1+(NW))=(2pr)*dl^2; %Z(i+1,j+1) end if i+1NW>0 Z(i,i+1(NW))=(2pr)*dl^2; %Z(i+1,j1) end if i 1*NW>0 Z(i,i1*(NW))=2*dw^2+2*(2pr)*dl^22*dw^2*dl; %Z(i,j1) end if i 2*NW>0 Z(i,i2*(NW))=dw^2; %Z(i,j2) end C(i)=Z(i,:)*[zeros((NL1)*NW+1,1);2*sin(3*pi*(i((NL1)*NW+1))/(NW1))*ones(NW2,1);0]; e nd %Spy on Z% spy(Z) %CALCULATION & REARRANGEMENT% z=Z\ C; for j=1:NL for i=1:NW Zsol(j,i)=z((j1)*NW+i,1); end end %PLOT OF SOLUTION% figure; x=0:dw:W; y=0:dl:L; [X,Y]=meshgrid(x,y); hold on %Upper Z Surface surf(X,Y,Zsol+H/2) %Lower Z Surface surf(X,Y,ZsolH/2) %Upper Y Surface [X1,Y1]=meshgrid(x,L); for i=1:NW Zsol1(i,:)=Zsol(NL,:)+H/2(H)*i/NW; end surf(X1,Y1,Zsol1) %Lower Y Surface [X2,Y2]=meshgrid(x,0); for i=1:NW Zsol2(i,:)=Zsol(1,:)+H/2(H)*i/NW; end PAGE 73 66 Appendix G Continued surf(X2,Y2,Zsol2) %Upper X Surface [X3,Y3]=meshgrid(W,y); for i=1:NL Zsol3(:,i)=Zsol(:,NW)+H/2(H)*i/(NL); end surf(X3,Y3,Zsol3) %Lower X Surface [X4,Y4]=meshgrid(0,y); for i=1:NL Zsol4(:,i)=Zsol(:,1)+H/2(H)*i/(NL); end surf(X4,Y4,Zsol4) colormap cool colorbar shading interp axis equal axis off zoom(1.4) lightangle(30,30) camorbit(155,72) set(gcf,'Renderer', 'zbuffer') set(findobj(gca,'type', 'surface'),... 'FaceLighting', 'phong', ... 'AmbientStrength',0.3,'DiffuseStrength',.8,... 'SpecularStrength',.9,'SpecularExponent',40,... 'BackFaceLighting', 'unlit') Determinant Behavior w.r.t. Nx This m file was used to determine how changing certain function values would give rise to a difference in the functions determinant. It was also used as a means to determine an appropriate initial guess for the critical magnitude of compression. function AnalyticalBuckling clc clear all close all m1=1; %Wavenumber W=30; %Width of Gel L=5; %Length of Gel Nx0=0; for i=1:85 Nx(i,1)=Nx0+10*(i1); D=10; %Flexural Rigidity alpham=m1*pi/W; PAGE 74 67 Appendix G Continued nu=0.25; %Poisson Ratio alpha=sqrt(m1^2*pi^2/(W^2)+sqrt(Nx(i,1)/D*m1^2*pi^2/(W^2))); beta=sqrt(m1^2*pi^2/(W^2)+sqrt(Nx(i,1)/D*m1^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y(i,1)=det(A); end m2=2; for i=1:46 Nx1(i,1)=Nx0+5*(i1); alpham=m2*pi/W; alpha=sqrt(m2^2*pi^2/(W^2)+sqrt(Nx1(i,1)/D*m2^2*pi^2/(W^2))); beta=sqrt(m2^2*pi^2/(W^2)+sqrt(Nx1(i,1)/D*m2^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y1(i,1)=det(A); end m3=3; for i=1:46 Nx2(i,1)=Nx0+2.5*(i1); alpham=m3*pi/W; alpha=sqrt(m3^2*pi^2/(W^2)+sqrt(Nx2(i,1)/D*m3^2*pi^2/(W^2))); beta=sqrt(m3^2*pi^2/(W^2)+sqrt(Nx2(i,1)/D*m3^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y2(i,1)=det(A); end PAGE 75 68 Appendix G Continued %Determinant Plot plot(Nx,y,'*' ,Nx1,y1,'d' ,Nx2,y2,'p' ) xlabel('N_x', 'Fontsize',14) ylabel('Determinant', 'Fontsize',14) title('Plot of Determinant against N_x', 'Fontsize',16) legend('m = 1', 'm = 2', 'm = 3') Determinant Behavior w.r.t. D This m file changes the value of the flexural rigidity and calculates the value of the determinant. This was used for initial guess purposes as well as a means to assess how changing the wavenumber changes the determinant behavior. function AnalyticalBuckling3 %Change D and find the Determinant clc clear all close all m1=1; %Wavenumber W=30; %Width of Gel L=10; %Length of Gel Nx=50; D0=5; %Flexural Rigidity for i=1:79 D1(i,1)=D0+2.5*(i1); alpham=m1*pi/W; nu=0.25; %Poisson Ratio alpha=sqrt(m1^2*pi^2/(W^2)+sqrt(Nx/D1(i,1)*m1^2*pi^2/(W^2))); beta=sqrt(m1^2*pi^2/(W^2)+sqrt(Nx/D1(i,1)*m1^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y1(i,1)=det(A); end D02=20; m2=2; PAGE 76 69 Appendix G Continued for i=1:73 D2(i,1)=D02+2.5*(i1); alpham=m2*pi/W; nu=0.25; %Poisson Ratio alpha=sqrt(m2^2*pi^2/(W^2)+sqrt(Nx/D2(i,1)*m2^2*pi^2/(W^2))); beta=sqrt(m2^2*pi^2/(W^2)+sqrt(Nx/D2(i,1)*m2^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y2(i,1)=det(A); end D03=40; m3=3; for i=1:65 D3(i,1)=D03+2.5*(i1); alpham=m3*pi/W; nu=0.25; %Poisson Ratio alpha=sqrt(m3^2*pi^2/(W^2)+sqrt(Nx/D3(i,1)*m3^2*pi^2/(W^2))); beta=sqrt(m3^2*pi^2/(W^2)+sqrt(Nx/D3(i,1)*m3^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; y3(i,1)=det(A); end %Determinant Plot plot(D1,y1,'*' ,D2,y2,'*' ,D3,y3,'*' ) xlabel('(D) Flexural Rigidity', 'Fontsize',14) ylabel('Determinant', 'Fontsize',14) title('Plot of Determinant against D', 'Fontsize',16) legend('m = 1', 'm = 2', 'm = 3') axis([0 200 13 8]) PAGE 77 70 Appendix G Continued Determinant Behavior w.r.t. L/W This m file analyzes how the length to width ratio affects the determinant. function AnalyticalBuckling6 clc clear all close all m=1; %Wavenumber W=30; %Width of Gel L0=3; Nx=20; nu=0.25; for i=1:38 L(i,1)=L0+.5*(i1); D=10; %Flexural Rigidity alpham=m*pi/W; alpha=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/D*m^2*pi^2/(W^2))); beta=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/D*m^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L(i,1))alpham^2*nu*exp(alpha*L(i,1))) ... (alpha^2*exp(alpha*L(i,1))alpham^2*nu*exp(alpha*L(i,1))) ... ( beta^2*cos(beta*L(i,1))alpham^2*nu*cos(beta*L(i,1))) ... ( beta^2*sin(beta*L(i,1))alpham^2*nu*sin(beta*L(i,1)));... ( alpha^3*exp(alpha*L(i,1))+alpham^2*(2nu)*alpha*exp(alpha*L(i,1))) ... (alpha^3*exp(alpha*L(i,1))alpham^2*(2nu)*alpha*exp(alpha*L(i,1))) ... (beta^3*sin(beta*L(i,1))+alpham^2*(2nu)*beta*sin(beta*L(i,1))) ... ( beta^3*cos(beta*L(i,1))alpham^2*(2nu)*beta*cos(beta*L(i,1)))]; y(i,1)=det(A); end %Determinant Plot plot(L/W,y,'*' ) xlabel('L/W', 'Fontsize',14) ylabel('Determinant', 'Fontsize',14) title('Plot of Determinant against L/W', 'Fontsize',16) PAGE 78 71 Appendix G Continued Buckling Solution and Critical Buckling Magnitude This m file solves for the critical buckling solution at different values of the wavenumber. This is done by solving the determinant expression for Nx and using the values that give rise to this solution to generate a plot of the critical buckling solution. function AnalyticalBuckling2 clc clear all close all %Generate Buckling Plot m=5; L=5; %Length in y direction W=30; %Width in x direction D=10; %Flexural Rigidity nu=0.25; %Poisson Ratio params.D=D; params.l=L; params.w=W; params.nu=nu; %Change the value of m and find smallest value of Nx for solution for i=1:37 mplot(i,1)=1+(i1)*0.25; if i==1 Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),50); else Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),Nxplot(i1,1)); end Nxplot(i,1)=Nxsol; if i>=2 if Nxplot(i,1)< Nxplot(i1,1) Nxmin=Nxplot(i,1); mmin=mplot(i,1); end end end disp(' ') disp(['Critical Nx = num2str(Nxmin)]) disp(' ') disp(['Critical Wavenumber = num2str(mmin)]) disp(' ') %Stability Plot plot(mplot,Nxplot,'*' ) xlabel('m', 'Fontsize',14) ylabel('N_x', 'Fontsize',14) PAGE 79 72 Appendix G Continued title('Plot of N_x against m', 'Fontsize',16) %Retrieve Constants from Constant Function Cval=buckling2(mmin,Nxmin,params); alpha=sqrt(mmin^2*pi^2/(W^2)+sqrt(Nxplot(17,1)/D*mmin^2*pi^2/(W^2))); beta=sqrt(mmin^2*pi^2/(W^2)+sqrt(Nxplot(17,1)/D*mmin^2*pi^2/(W^2))); [X,Y]=meshgrid(0:0.75:W,0:0.75:L); Z=sin(mmin.*pi.*X./W).*(Cval(1,1).*exp(alpha.*Y)+Cval(2,1).*exp(alpha.*Y)+Cval(3,1).*cos(beta.*Y)+Cval(4,1).*s in(beta.*Y)); %Plot of Buckling Solution figure; surf(X,Y,Z) colormap cool colorbar shading interp axis equal axis off lightangle(30,30) camorbit(105,3) set(gcf,'Renderer', 'zbuffer') set(findobj(gca,'type', 'surface'),... 'FaceLighting', 'phong', ... 'AmbientStrength',0.3,'DiffuseStrength',.8,... 'SpecularStrength',.9,'SpecularExponent',40,... 'BackFaceLighting', 'unlit') %Determinant Function function f=buckling(m,Nx,params) l=params.l; %Length of Gel w=params.w; %Width of gel alpham=m*pi/w; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A=1; B=1; C=1; D=0; E=alpha; F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); PAGE 80 73 Appendix G Continued %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); %Constant Function function C=buckling2(m,Nx,params) L=params.l; %Length of Gel W=params.w; %Width of gel alpham=m*pi/W; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/Def*m^2*pi^2/(W^2))); beta=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/Def*m^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; c=[0;0;0;0]; Amat=[A c]; %Making C4=1 Aref=rref(Amat); C=[Aref(1,4);Aref(2,4);Aref(3,4);1]; Change m and L/W to Determine Critical Nx function AnalyticalBuckling9 clc clear all close all %Generate Buckling Plot L=7.5; %Length in y direction W=30; %Width in x direction D=10; %Flexural Rigidity nu=0.25; %Poisson Ratio params.D=D; params.l=L; params.w=W; params.nu=nu; %Change the value of m and find smallest value of Nx for solution for i=1:37 PAGE 81 74 Appendix G Continued mplot(i,1)=1+(i1)*0.25; if i==1 Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),80); else Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),Nxplot(i1,1)); end Nxplot(i,1)=Nxsol; end L1=15; params.l=L1; for i=1:37 mplot1(i,1)=1+(i1)*0.25; if i==1 Nxsol1=fzero(@(Nx)buckling(mplot1(i,1),Nx,params),50); else Nxsol1=fzero(@(Nx)buckling(mplot1(i,1),Nx,params),Nxplot1(i1,1)); end Nxplot1(i,1)=Nxsol1; end L2=30; params.l=L2; for i=1:37 mplot2(i,1)=1+(i1)*0.25; if i==1 Nxsol2=fzero(@(Nx)buckling(mplot2(i,1),Nx,params),50); else Nxsol2=fzero(@(Nx)buckling(mplot2(i,1),Nx,params),Nxplot2(i1,1)); end Nxplot2(i,1)=Nxsol2; end %Stability Plot plot(mplot,Nxplot,'*' ,mplot1,Nxplot1,'d' ,mplot2,Nxplot2,'p' ) xlabel('Wavenumber (m)', 'Fontsize',14) ylabel('N_x', 'Fontsize',14) title('Plot of N_x against m', 'Fontsize',16) legend('L/W = 0.25', 'L/W = 0.5', 'L/W = 1') %Determinant Function function f=buckling(m,Nx,params) l=params.l; %Length of Gel w=params.w; %Width of gel alpham=m*pi/w; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A=1; B=1; C=1; D=0; E=alpha; PAGE 82 75 Appendix G Continued F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); Change m and D to Determine Critical Nx function AnalyticalBuckling8 clc clear all close all %Generate Buckling Plot m=5; L=5; %Length in y direction W=30; %Width in x direction D=10; %Flexural Rigidity nu=0.25; %Poisson Ratio params.D=D; params.l=L; params.w=W; params.nu=nu; %Change the value of m and find smallest value of Nx for solution for i=1:37 mplot(i,1)=1+(i1)*0.25; if i==1 Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),50); else Nxsol=fzero(@(Nx)buckling(mplot(i,1),Nx,params),Nxplot(i1,1)); end Nxplot(i,1)=Nxsol; end D1=20; params.D=D1; for i=1:37 mplot1(i,1)=1+(i1)*0.25; if i==1 PAGE 83 76 Appendix G Continued Nxsol1=fzero(@(Nx)buckling(mplot1(i,1),Nx,params),50); else Nxsol1=fzero(@(Nx)buckling(mplot1(i,1),Nx,params),Nxplot1(i1,1)); end Nxplot1(i,1)=Nxsol1; end D2=30; params.D=D2; for i=1:37 mplot2(i,1)=1+(i1)*0.25; if i==1 Nxsol2=fzero(@(Nx)buckling(mplot2(i,1),Nx,params),50); else Nxsol2=fzero(@(Nx)buckling(mplot2(i,1),Nx,params),Nxplot2(i1,1)); end Nxplot2(i,1)=Nxsol2; end %Stability Plot plot(mplot,Nxplot,'*' ,mplot1,Nxplot1,'d' ,mplot2,Nxplot2,'p' ) xlabel('Wavenumber (m)', 'Fontsize',14) ylabel('N_x', 'Fontsize',14) title('Plot of N_x against m', 'Fontsize',16) legend('D = 10', 'D = 20', 'D = 30') %Determinant Function function f=buckling(m,Nx,params) l=params.l; %Length of Gel w =params.w; %Width of gel alpham=m*pi/w; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A=1; B=1; C=1; D=0; E=alpha; F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... PAGE 84 77 Appendix G Continued B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); Change m and Nx to Determine C ritical D This m file calculates the minimum critical value of flexural rigidity that will give rise to a buckling solution at different wavenumbers. function AnalyticalBuckling4 clc clear all close all L=5; %Length in y direction W=30; %Width in x direction Nx=30; nu=0.25; %Poisson Ratio params.Nx=Nx; params.l=L; params.w=W; params.nu=nu; %Change the value of m and find smallest value of D for solution for i=1:37 mplot(i,1)=1+(i1)*0.25; if i==1 Dsol=fzero(@(Def)buckling(mplot(i,1),Def,params),20); else Dsol=fzero(@(Def)buckling(mplot(i,1),Def,params),Dplot(i1,1)); end Dplot(i,1)=Dsol; end Nx=40; params.Nx=Nx; for i=1:37 mplot1(i,1)=1+(i1)*0.25; if i==1 Dsol1=fzero(@(Def)buckling(mplot1(i,1),Def,params),20); else Dsol1=fzero(@(Def)buckling(mplot1(i,1),Def,params),Dplot1(i1,1)); end Dplot1(i,1)=Dsol1; end Nx=50; params.Nx=Nx; for i=1:37 mplot2(i,1)=1+(i1)*0.25; PAGE 85 78 Appendix G Continued if i==1 Dsol2=fzero(@(Def)buckling(mplot2(i,1),Def,params),20); else Dsol2=fzero(@(Def)buckling(mplot2(i,1),Def,params),Dplot2(i1,1)); end Dplot2(i,1)=Dsol2; end %Stability Plot plot(mplot,Dplot,'*' ,mplot1,Dplot1,'d' ,mplot2,Dplot2,'p' ) xlabel('Wavenumber (m)', 'Fontsize',14) ylabel('Flexural Rigidity (D)', 'Fontsize',14) title('Plot of D against m', 'Fontsize',16) legend('N_x= 30', 'N_x= 40', 'N_x= 50') %Determinant Function function f=buckling(m,Def,params) l=params.l; %Length of Gel w=params.w; %Width of gel alpham=m*pi/w; nu=params.nu; %Poisson Ratio Nx=params.Nx; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A =1; B=1; C=1; D=0; E=alpha; F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); PAGE 86 79 Appendix G Continued Change m and Nx to Determine C riti cal L/W This m file calculates the minimum critical value of the length and width ratio that will give rise to a buckling solution at different wavenumbers. function AnalyticalBuckling7 clc clear all close all W=30; %Width in x direction Def=10; %Flexural Rigidity nu=0.25; %Poisson Ratio Nx=10; params.Nx=Nx; params.nu=nu; params.w=W; params.Def=Def; %Change the value of m and find smallest value of L for solution for i=1:33 mplot(i,1)=1+(i1)*0.25; if i==1 Lsol=fzero(@(L)buckling(mplot(i,1),L,params),2); else Lsol=fzero(@(L)buckling(mplot(i,1),L,params),Lplot(i1,1)); end Lplot(i,1)=Lsol; end Nx=20; params.Nx=Nx; for i=1:48 mplot1(i,1)=1+(i1)*0.25; if i==1 Lsol1=fzero(@(L)buckling(mplot1(i,1),L,params),2); else Lsol1=fzero(@(L)buckling(mplot1(i,1),L,params),Lplot1(i1,1)); end Lplot1(i,1)=Lsol1; end Nx=30; params.Nx=Nx; for i=1:30 mplot2(i,1)=1+(i1)*0.5; if i==1 Lsol2=fzero(@(L)buckling(mplot2(i,1),L,params),2); else Lsol2=fzero(@(L)buckling(mplot2(i,1),L,params),Lplot2(i1,1)); PAGE 87 80 Appendix G Continued end Lplot2(i,1)=Lsol2; end %Stability Plot plot(mplot,Lplot/W,'*' ,mplot1,Lplot1/W,'d' ,mplot2,Lplot2/W,'p' ) xlabel('m', 'Fontsize',14) ylabel('L / W', 'Fontsize',14) title('Plot of L / W against m', 'Fontsize',16) legend('N_x= 10', 'N_x= 20', 'N_x= 30') %Determinant Function function f=buckling(m,l,params) nu=params.nu; %Poisson Ratio w=params.w; %Width of gel alpham=m*pi/w; Def=params.Def; %Flexural Rigidity Nx=params.Nx; alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A=1; B=1; C=1; D=0; E=alpha; F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); Strain Energy and Work Equation Generator This m file uses the Symbolic Math Toolbox in MATLAB to determine expressions for the work and strain energy of the hydrogel. PAGE 88 81 Appendix G Continued clc clear all syms c1 c2 c3 c4 a b x y alpham alpha beta nu D Nx W L %Solution to Biharmonic Equation Z=sin(alpham*x)*(c1*exp(alpha*y)+c2*exp(alpha*y)+c3*cos(beta*y)+c4*sin(beta*y)) %Equation for Work Done on Hydrogel Work=Nx/2*int(int((diff(Z,x)^2),x,0,W),y,0,L); %Simplified Work Expression SimpleWork=simplify(Work) %Differentials in Strain Energy Equation difftwox=diff(diff(Z,x),x); difftwoy=diff(diff(Z,y),y); difftwoxy=diff(diff(Z,y),x); %Strain Energy Equation StrainEnergy=D/2*int(int((difftwox+difftwoy)^22*(1nu)*(difftwox*difftwoy... (difftwoxy)^2),y,0,L),x,0,W); %Simplified Strain Energy Expression SimpleStrainEnergy=simplify(StrainEnergy) Strain Energy and Work Function This m file uses the determined expressions for work and strain energy to determine critical values of Nx at different wavenumbers. This is used as a means to ver ify that the results obtained are in agreement even though different methods were required to obtain them. function StrainEnergy clc close all clear all W=30; L=5; D=10; nu=0.25; params.W=W; params.L=L; PAGE 89 82 Appendix G Continued params.nu=nu; params.D=D; N=18; for i=1:N mval(i,1)=1+(i1)*0.5; if i==1 Nxsol(i,1)=fzero(@(Nx)buckling(mval(i,1),Nx,params),30); else Nxsol(i,1)=fzero(@(Nx)buckling(mval(i,1),Nx,params),Nxsol(i1,1)); end end %Value of m to be analyzed meval=1; %Constants that cause buckling C=buckling2(meval,Nxsol(meval+(meval1),1),params); params.c1=C(1); params.c2=C(2); params.c3=C(3); params.c4=C(4); for i=1:N+10 Nxval(i,1)=0.5+(i1)*1.25; work1(i,1)=workfun(meval,Nxval(i,1),params); strain1(i,1)=strainfun(meval,Nxval(i,1),params); end plot(Nxval,work1,'*' ,Nxval,strain1,'^r') xlabel('N_x', 'Fontsize',14) ylabel('Energy', 'Fontsize',14) title(['Plot of Energy against N_x when m = num2str(meval)],'Fontsize',16) legend('Work', 'Strain') Nxcrit=fzero(@(Nx)(workfun(meval,Nx,params)(strainfun(meval,Nx,params))),30) function f=workfun(m,Nx,params) c1=params.c1; c2=params.c2; c3=params.c3; c4=params.c4; D=params.D; W=params.W; L=params.L; alpham=m*pi/W; alpha=sqrt(alpham^2+sqrt(Nx/D*alpham^2)); beta=sqrt(alpham^2+sqrt(Nx/D*alpham^2)); f=1/64*Nx*alpham*(4*c2^2*sin(2*alpham*W)*beta*alpha^2+8*c1*c4*... exp(alpha*L)*alpha*beta^2*sin(beta*L2*alpham*W)8*c1^2*alpham*W*... PAGE 90 83 Appendix G Continued exp(alpha*L)^2*beta*alpha^28*c2^2*alpham*W*beta^3+8*c2^2*alpham*... W*exp(alpha*L)^2*beta^3+8*c3^2*alpham*W*L*alpha^3*beta+2*c3*c4*... sin(2*beta*L2*alpham*W)*alpha*beta^28*c1*c3*exp(alpha*L)*alpha^2*... beta*sin(beta*L+2*alpham*W)+8*c3*c4*alpham*W*alpha*beta^2+8*c3*c4*... alpham*W*alpha^38*c2^2*alpham*W*beta*alpha^2+8*c1^2*alpham*W*beta*... alpha^2+4*c1^2*sin(2*alpham*W)*beta*alpha^2+8*c1^2*alpham*W*beta^3+4*.. c3*c4*sin(2*alpham*W)*alpha^34*c2^2*sin(2*alpham*W)*beta^3+4*c1^2*... sin(2*alpham*W)*beta^3+16*c1*c3*alpha^2*beta*sin(2*alpham*W)+4*c3*c4*.. sin(2*alpham*W)*alpha*beta^2+16*c1*c4*alpha*beta^2*sin(2*alpham*W)... 32*c2*c3*alpham*W*alpha^2*beta+32*c1*c4*alpham*W*alpha*beta^2+32*c1*... c3*alpham*W*alpha^2*beta+32*c2*c4*alpham*W*alpha*beta^2+16*c2*c4*... alpha*beta^2*sin(2*alpham*W)16*c2*c3*alpha^2*beta*sin(2*alpham*W)... c3^2*cos(2*beta*L+2*alpham*W)*alpha^3c4^2*cos(2*beta*L2*alpham*W)... *alpha^3+c3^2*cos(2*beta*L2*alpham*W)*alpha^3+c4^2*cos(2*beta*L+2*... alpham*W)*alpha^3+8*c3^2*alpham*W*L*alpha*beta^3+8*c4^2*alpham*W*L*... alpha^3*beta+8*c4^2*alpham*W*L*alpha*beta^3+32*c1*c2*alpham*W*L*... alpha^3*beta+32*c1*c2*alpham*W*L*alpha*beta^34*c1^2*sin(2*alpham*W)... *exp(alpha*L)^2*beta^3c4^2*cos(2*beta*L2*alpham*W)*alpha*... beta^2+c3^2*cos(2*beta*L2*alpham*W)*alpha*beta^2+4*c2^2*... sin(2*alpham*W)*exp(alpha*L)^2*beta^32*c3*c4*sin(2*beta*L+2*... alpham*W)*alpha^3c3^2*cos(2*beta*L+2*alpham*W)*alpha*beta^2+c4^2*... cos(2*beta*L+2*alpham*W)*alpha*beta^2+2*c3*c4*sin(2*beta*L2*alpham*... W)*alpha^3+4*c3^2*alpham*W*sin(2*beta*L)*alpha^38*c2*c3*... exp(alpha*L)*alpha*beta^2*cos(beta*L+2*alpham*W)+8*c2^2*alpham*W*... exp(alpha*L)^2*beta*alpha^2+16*c1*c2*sin(2*alpham*W)*L*alpha^3*beta... 8*c1*c3*exp(alpha*L)*alpha*beta^2*cos(beta*L+2*alpham*W)+4*c3^2*... alpham*W*sin(2*beta*L)*alpha*beta^24*c1^2*sin(2*alpham*W)*... exp(alpha*L)^2*beta*alpha^2+4*c2^2*sin(2*alpham*W)*exp(alpha*L)^2*... beta*alpha^22*c3*c4*sin(2*beta*L+2*alpham*W)*alpha*beta^28*c1^2*... alpham*W*exp(alpha*L)^2*beta^34*c4^2*alpham*W*sin(2*beta*L)*... alpha^34*c4^2*alpham*W*sin(2*beta*L)*alpha*beta^2+16*c1*c2*... PAGE 91 84 Appendix G Continued sin(2*alpham*W)*L*alpha*beta^3+4*c4^2*sin(2*alpham*W)*L*alpha^3*... beta+4*c4^2*sin(2*alpham*W)*L*alpha*beta^3+4*c3^2*sin(2*alpham*W)*L*... alpha^3*beta+4*c3^2*sin(2*alpham*W)*L*alpha*beta^38*c3*c4*alpham*W*... cos(2*beta*L)*alpha^3 8*c3*c4*alpham*W*cos(2*beta*L)*alpha*beta^28*... c2*c4*exp(alpha*L)*alpha^2*beta*cos(beta*L+2*alpham*W)8*c2*c4*... exp(alpha*L)*alpha*beta^2*sin(beta*L+2*alpham*W)8*c1*c4*exp(alpha*L)... *alpha^2*beta*cos(beta*L2*alpham*W)+8*c1*c3*exp(a lpha*L)*alpha*... beta^2*cos(beta*L2*alpham*W)+8*c1*c3*exp(alpha*L)*alpha^2*beta... *sin(beta*L2*alpham*W)+8*c1*c4*exp(alpha*L)*alpha^2*beta*... cos(beta*L+2*alpham*W)8*c1*c4*exp(alpha*L)*alpha*beta^2*... sin(beta*L+2*alpham*W)+8*c2*c3*exp(alpha*L)*alpha^2*beta*... sin(beta*L+2*alpham*W)+8*c2*c3*exp(alpha*L)*alpha*beta^2*... cos(beta*L2*alpham*W)8*c2*c3*exp(alpha*L)*alpha^2*beta*sin(beta*L... 2*alpham*W)+8*c2*c4*exp(alpha*L)*alpha^2*beta*cos(beta*L2*alpham*W)... +8*c2*c4*exp(alpha*L)*alpha*beta^2*sin(beta*L2*alpham*W)32*c1*c4*... alpham*W*exp(alpha*L)*alpha*beta^2*cos(beta*L)32*c1*c4*alpham*W*... exp(alpha*L)*alpha^2*beta*sin(beta*L)32*c1*c3*alpham*W*... exp(alpha*L)*alpha^2*beta*cos(beta*L)+32*c1*c3*alpham*W*... exp(alpha*L)*alpha*beta^2*sin(beta*L)+32*c2*c3*alpham*W*... exp(alpha*L)*alpha^2*beta*cos(beta*L)+32*c2*c3*alpham*W*exp(alpha*L)*.. alpha*beta^2*sin(beta*L)32*c2*c4*alpham*W*exp(alpha*L)*alpha*beta^2*... cos(beta*L)+32*c2*c4*alpham*W*exp(alpha*L)*alpha^2*beta*sin(beta*L))... /(alpha^2+beta^2)/alpha/beta; function f=strainfun(m,Nx,params) D=params.D; c1=params.c1; c2=params.c2; c3=params.c3; c4=params.c4; W=params.W; L=params.L; nu=params.nu; alpham=m*pi/W; alpha=sqrt(alpham^2+sqrt(Nx/D*alpham^2)); beta=sqrt(alpham^2+sqrt(Nx/D*alpham^2)); f=1/64*D*(32*alpham^5*W*c2*c4*beta^2*alpha+16*alpha*L*alpham^3*c4^2*... beta^5*W+16*alpha*L*alpham^3*c3^2*beta^5*W+8*alpha*L*alpham^5*W*c4^2*.. PAGE 92 85 Appendix G Continued beta^3+8*alpha*L*W*c4^2*beta^7*alpham8*alpham^2*alpha^4*beta*nu*c3*... c2*exp(alpha*L)*sin(beta*L2*alpham*W)+8*alpham^2*alpha^4*beta*nu*c3*... c1*exp(alpha*L)*sin(beta*L2*alpham*W)8*alpham^4*c1*c4*... exp(alpha*L)*alpha*beta^2*sin(beta*L2*alpham*W)16*alpham^2*... alpha^2*beta^3*c2*c4*exp(alpha*L)*cos(beta*L+2*alpham*W)8*c1^2*... alpha^6*exp(alpha*L)^2*beta*W*alpham32*c2*alpha^3*c3*beta^4*... exp(alpha*L)*sin(beta*L)*W*alpham+8*beta^4*W*c3*c4*alpham*alpha^316*... beta^2*alpham^3*W*c3*c4*alpha^3+32*beta^2*alpham^3*nu*c3*c4*W*... alpha^3+8*alpham^5*W*c3*c4*alpha^3+8*alpha*L*alpham^5*W*c3^2*beta^3+... 32*alpha*L*alpham^5*c1*c2*W*beta^3+32*beta^4*alpham^3*nu*c3*c4*W*... alpha+8*beta^6*W*c3*c4*alpham*alpha8*alpha^6*c2^2*W*beta*alpham16*... alpham^4*sin(2*alpham*W)*beta^2*c1*c4*alpha+8*c1^2*alpha^6*W*beta*... alpham4*beta^4*c4^2*sin(2*beta*L)*alpha^3*W*alpham+32*c2*alpha^3*c4*... beta^4*exp(alpha*L)*cos(beta*L)*W*alpham24*alpham^2*alpha^2*beta^3*... nu*c1*c3*exp(alpha*L)*sin(beta*L2*alpham*W)32*alpham^5*c1*c4*... exp(alpha*L)*alpha^2*beta*sin(beta*L)*W32*alpham^5*c1*c4*... exp(alpha*L)*alpha*beta^2*cos(beta*L)*W8*alpham^4*c2*c4*... exp(alpha*L)*alpha^2*beta*cos(beta*L2*alpham*W)+8*alpha*L*W*c3^2*... beta^7*alpham32*alpham^3*c1*alpha^4*nu*W*c3*beta+16*sin(2*alpham*W)*... c1*alpha^4*c3*beta^3+16*alpham^3*W*c1^2*alpha^2*beta^38*alpham^5* ... c2^2*W*beta*alpha^24*alpham^4*c2^2*exp(alpha*L)^2*beta*alpha^2*... sin(2*alpham*W)+8*alpham^2*alpha*beta^4*nu*c1*c3*exp(alpha*L)*... cos(beta*L+2*alpham*W)8*alpham^2*alpha^4*beta*nu*c4*c2*exp(alpha*L)*... cos(beta*L+2*alpham*W)32*alpham^3*alpha^3*beta^2*nu*c1*c4*... exp(alpha*L)*cos(beta*L)*W+8*alpham^5*c2^2*exp(alpha*L)^2*beta^3*W+... 32*alpham^5*c1*alpha^2*W*c3*beta+4*alpha^6*sin(2*alpham*W)*c2^2*... beta+4*alpham^4*sin(2*alpham*W)*c2^2*beta^3+32*alpham^3*nu*c2*c3*... alpha^4*W*beta+32*alpha^4*alpham^3*c2^2*W*nu*beta+32*alpha^2*alpham^... 3*c2^2*W*nu*beta^34*c1^2*alpha^4*sin(2*alpham*W)*beta^34*c4^2*... alpham^5*sin(2*beta*L)*alpha*beta^2*W+32*alpham^3*alpha^4*beta*nu*... c3*c1*exp(alpha*L)*cos(beta*L)*W32*alpham^3*alpha^4*beta*nu*c3*c2*... exp(alpha*L)*cos(beta*L)*W8*alpham^4*c2*c4*exp(alpha*L)*alpha*beta^... 2*sin(beta*L2*alpham*W)+16*alpham^2*alpha^3*beta^2*c3*c2*... exp(alpha*L)*cos(beta*L+2*alpham*W)+8*alpham^2*alpha^4*beta*nu*c4*... c1*exp(alpha*L)*cos(beta*L+2*alpham*W)16*alpham^3*alpha*beta^4*nu*... PAGE 93 86 Appendix G Continued c4^2*sin(2*beta*L)*W+8*alpham^5*W*c3*c4*alpha*beta^2+16*sin(2*alpham*.. W)*beta^4*c1*alpha^3*c416*sin(2*alpham*W)*c2*alpha^4*c3*beta^3+16*... alpham^2*sin(2*alpham*W)*c1*alpha^4*nu*c3*beta32*c1*alpha^4*beta^3*... W*c3*alpham+32*alpham^3*nu*c2*c3*beta^3*W*alpha^2+48*alpham^2*... sin(2*alpham*W)*nu*c2*alpha^3*c4*beta^232*alpham^2*sin(2*alpham*W)*... c2*alpha^3*c4*beta^248*alpham^2*sin(2*alpham*W)*beta^3*c1*nu*... alpha^2*c38*alpham^5*c2^2*W*beta^364*alpham^3*c2*alpha^2*c3*beta^3*... W+32*alpham^3*nu*c2*alpha^3*c4*beta^2*W8*alpha^2*alpham^2*... sin(2*alpham*W)*c2^2*beta^3+32*alpham^2*sin(2*alpham*W)*beta^3*c1*... alpha^2*c3+48*alpham^2*sin(2*alpham*W)*beta^2*c1*nu*alpha^3*c4+4*... c3^2*alpham^5*sin(2*beta*L)*alpha^3*W+24*alpham^2*alpha^3*beta^2*nu*... c1*c3*exp(alpha*L)*cos(beta*L2*alpham*W)+4*c3^2*alpham^5*... sin(2*beta*L)*alpha*beta^2*W8*alpham^5*c3*c4*cos(2*beta*L)*alpha*... beta^2*W+8*alpham^4*c2*c4*exp(alpha*L)*alpha^2*beta*... cos(beta*L+2*alpham*W)+16*alpham^2*alpha^2*beta^3*c1*c4*... exp(alpha*L)*cos(beta*L+2*alpham*W)+c3^2*alpham^4*alpha*beta^2*... cos(2*beta*L+2*alpham*W)+32*c1*alpha^3*c4*beta^4*exp(alpha*L)*... cos(beta*L)*W*alpham8*alpha^4*c2^2*W*beta^3*alpham4*alpha*L*... alpham^4*sin(2*alpham*W)*c3^2*beta^3+8*c2*alpha^3*c3*beta^4*... exp(alpha*L)*cos(beta*L2*alpham*W)2*c3*beta^4*c4*alpha^3*... sin(2*beta*L2*alpham*W)8*c2*alpha^3*c4*beta^4*exp(alpha*L)*... sin(beta*L+2*alpham*W)2*c3*beta^6*c4*alpha*sin(2*beta*L2*... alpham*W)+c3^2*alpham^4*alpha^3*cos(2*beta*L+2*alpham*W)32*c1^2*... alpha^4*alpham^3*W*nu*beta32*alpham^3*nu*c1^2*alpha^2*W*beta^3+4*... alpha^4*sin(2*alpham*W)*c2^2*beta^316*alpham^2*alpha^2*beta^3*c1*... c3*exp(alpha*L)*sin(beta*L+2*alpham*W)+4*alpham^2*c3*c4*beta^2*... alpha^3*sin(2*beta*L+2*alpham*W)+2*beta^2*c4^2*alpham^2*alpha^3*... cos(2*beta*L2*alpham*W)2*alpham^2*alpha*beta^4*c3^2*... cos(2*beta*L2*alpham*W)+24*alpham^2*alpha^2*beta^3*nu*c2*c3*... exp(alpha*L)*sin(beta*L2*alpham*W)+64*alpham^3*alpha^2*beta^3*c2*... c4*exp(alpha*L)*sin(beta*L)*W2*alpham^2*alpha^3*beta^2*c3^2*... cos(2*beta*L2*alpham*W)4*L*sin(2*alpham*W)*beta^7*c3^2*... alphabeta^4*c3^2*alpha^3*cos(2*beta*L2*alpham*W)32*alpham^3*alpha*... beta^4*nu*c1*c4*exp(alpha*L)*cos(beta*L)*W+24*alpham^2*alpha^2*... beta^3*nu*c2*c4*exp(alpha*L)*cos(beta*L+2*alpham*W)32*alpham^3*... alpha^2*beta^3*nu*c2^2*exp(alpha*L)^2*W16*alpham^3*alpha^2*beta^3*... c1^2*exp(alpha*L)^2*W+24*alpham^2*alpha^3*beta^2*nu*c3*c2*... exp(alpha*L)*cos(beta*L2*alpham*W)8*alpham^5*c3*c4*cos(2*beta*L)*... PAGE 94 87 Appendix G Continued alpha^3*W+8*alpham^5*W*c1^2*beta^332*alpham^5*c1*c3*exp(alpha*L)*... alpha^2*beta*cos(beta*L)*W+32*alpham^5*c1*c3*exp(alpha*L)*alpha*... beta^2*sin(beta*L)*W+4*alpham^2*c3*c4*beta^4*alpha*... sin(2*beta*L+2*alpham*W)8*alpham^2*alpha*beta^4*nu*c2*c4*... exp(alpha*L)*sin(beta*L2*alpham*W)8*alpham^4*c1*c3*exp(alpha*L)*... alpha^2*beta*sin(beta*L2*alpham*W)24*alpham^2*alpha^2*beta^3*nu*... c2*c4*exp(alpha*L)*cos(beta*L2*alpham*W)+8*alpham^4*c2*c3*... exp(alpha*L)*alpha*beta^2*cos(beta*L+2*alpham*W)+8*alpham^2*c1^2*... alpha^2*sin(2*alpham*W)*beta^332*alpham^2*sin(2*alpham*W)*c2*... alpha^2*c3*beta^332*alpham^3*beta^3*c1*nu*alpha^2*c3*W32*beta^4*... c1*W*alpha^3*c4*alpham64*alpham^3*c1*alpha^3*beta^2*W*c4+32*alpha^... 7*L*c1*c2*W*alpham*beta64*alpha^5*L*alpham^3*c1*c2*W*beta+16*alpham^... 4*sin(2*alpham*W)*c2*c3*alpha^2*beta+32*alpham^5*beta^2*c1*W*c4*... alpha+32*alpha^5*L*c1*c2*W*alpham*beta^316*beta^4*alpham^3*W*c3*... c4*alpha+32*beta^4*alpham^3*c1*nu*c4*W*alpha24*alpham^2*alpha^3*... beta^2*nu*c1*c3*exp(alpha*L)*cos(beta*L+2*alpham*W)+4*alpham^4*... c1^2*exp(alpha*L)^2*beta*alpha^2*sin(2*alpham*W)64*alpham^3*c1*c4*... beta^3*exp(alpha*L)*alpha^2*sin(beta*L)*W+8*c2*alpha^3*c4*beta^4*... exp(alpha*L)*sin(beta*L2*alpham*W)2*alpham^4*c3*c4*alpha^3*... sin(2*beta*L2*alpham*W)+8*alpham^2*alpha*beta^4*nu*c1*c4*... exp(alpha*L)*sin(beta*L+2*alpham*W)+8*alpham^5*W*c1^2*beta*alpha^2... 4*alpham^4*c1^2*sin(2*alpham*W)*beta^3+32*alpham^3*nu*c2*c4*beta^4*... W*alpha+48*alpham^2*sin(2*alpham*W)*nu*c3*c2*beta^3*alpha^2+4*alpham^.. 4*sin(2*alpham*W)*c2^2*beta*alpha^24*c2^2*alpha^4*exp(alpha*L)^2*... beta^3*sin(2*alpham*W)+4*c1^2*alpha^6*exp(alpha*L)^2*beta*... sin(2*alpham*W)32*W*c2*alpha^3*c4*beta^4*alpham32*alpham^2*... sin(2*alpham*W)*c1*alpha^3*c4*beta^28*alpha^4*alpham^2*... sin(2*alpham*W)*c2^2*beta16*alpham^2*sin(2*alpham*W)*nu*c3*c2*... alpha^4*beta8*alpham^2*alpha^4*beta*nu*c4*c1*exp(alpha*L)*... cos(beta*L2*alpham*W)c4^2*alpham^4*alpha*beta^2*... cos(2*beta*L+2*alpham*W)+32*alpham^3*alpha^4*beta*nu*c4*c1*... exp(alpha*L)*sin(beta*L)*Wc3^2*alpham^4*alpha*beta^2*... cos(2*beta*L2*alpham*W)32*alpham^3*alpha*beta^4*nu*c3*c4*... cos(2*beta*L)*W+64*alpham^3*alpha^3*beta^2*c1*c4*exp(alpha*L)*... cos(beta*L)*W8*alpham^2*alpha*beta^4*nu*c1*c3*exp(alpha*L)*... cos(beta*L2*alpham*W)+8*c1^2*alpha^4*alpham^2*sin(2*alpham*W)*beta... PAGE 95 88 Appendix G Continued 2*beta^4*c4^2*alpham^2*alpha*cos(2*beta*L+2*alpham*W)+24*alpham^2*... alpha^2*beta^3*nu*c1*c4*exp(alpha*L)*cos(beta*L2*alpham*W)+64*... alpham^3*alpha^3*beta^2*c4*c2*exp(alpha*L)*cos(beta*L)*W+16*alpham^2*.. alpha^3*beta^2*c3*c1*exp(alpha*L)*cos(beta*L+2*alpham*W)8*c2*alpha^... 4*c4*beta^3*exp(alpha*L)*cos(beta*L+2*alpham*W)24*alpham^2*alpha^... 3*beta^2*nu*c3*c2*exp(alpha*L)*cos(beta*L+2*alpham*W)24*alpham^2*... alpha^3*beta^2*nu*c1*c4*exp(alpha*L)*sin(beta*L+2*alpham*W)+beta^4*... c3^2*alpha^3*cos(2*beta*L+2*alpham*W)+8*c1^2*alpha^4*W*beta^3*... alpham16*alpham^2*alpha^2*beta^3*c1*c4*exp(alpha*L)*... cos(beta*L2*alpham*W)+16*alpham^3*c2^2*alpha^4*exp(alpha*L)^2*... beta*W32*alpham^3*alpha^3*beta^2*nu*c3*c4*cos(2*beta*L)*W2*beta^2*... c4^2*alpham^2*alpha^3*cos(2*beta*L+2*alpham*W)+2*alpham^4*c3*c4*... alpha^3*sin(2*beta*L+2*alpham*W)+16*alpham^3*alpha^3*beta^2*c3*c4*... cos(2*beta*L)*W+64*alpham^3*alpha^2*beta^3*c2*c3*exp(alpha*L)*... cos(beta*L)*W16*alpham^4*sin(2*alpham*W)*c2*c4*beta^2*alpha+32*... beta^2*alpham^3*c1*nu*c4*alpha^3*W+16*alpha^3*L*alpham^3*c3^2*beta^3*.. W+64*alpham^3*beta^3*c1*alpha^2*c3*W16*alpham^2*sin(2*alpham*W)*nu*... c2*c4*beta^4*alpha+16*alpha^3*L*alpham^3*c4^2*beta^3*W16*L*... sin(2*alpham*W)*alpha^7*c2*c1*beta16*L*sin(2*alpham*W)*alpha^5*c2*... c1*beta^3+64*alpha^5*L*alpham^2*sin(2*alpham*W)*nu*c2*c1*beta32*... alpha^5*L*alpham^2*sin(2*alpham*W)*c2*c1*beta8*c1*alpha^3*c3*beta^4*... exp(alpha*L)*cos(beta*L+2*alpham*W)8*c2*alpha^4*c3*beta^3*... exp(alpha*L)*sin(beta*L2*alpham*W)4*c3*beta^4*c4*alpha^3*... sin(2*alpham*W)24*alpham^2*alpha^3*beta^2*nu*c4*c2*exp(alpha*L)*... sin(beta*L+2*alpham*W)+16*alpham^2*alpha^3*beta^2*c4*c2*exp(alpha*L)*.. sin(beta*L+2*alpham*W)+24*alpham^2*alpha^3*beta^2*nu*c4*c2*... exp(alpha*L)*sin(beta*L2*alpham*W)8*alpham^5*c1^2*... exp(alpha*L)^2*beta*alpha^2*W+4*alpham^4*c1^2*exp(alpha*L)^2*... beta^3*sin(2*alpham*W)8*alpham^2*alpha^4*beta*nu*c3*c1*... exp(alpha*L)*sin(beta*L+2*alpham*W)64*alpham^3*c2*alpha^3*c4*... beta^2*W+8*alpha^3*L*alpham^5*W*c4^2*beta+8*alpha^3*L*W*c4^2*beta^5*... alpham+32*alpham^5*c2*c3*exp(alpha*L)*alpha^2*beta*cos(beta*L)*W+24*... alpham^2*alpha^2*beta^3*nu*c1*c3*exp(alpha*L)*... sin(beta*L+2*alpham*W)+16*alpham^2*alpha^3*beta^2*c1*c4*exp(alpha*L)*... PAGE 96 89 Appendix G Continued sin(beta*L+2*alpham*W)16*alpham^2*alpha^2*beta^3*c2*c3*exp(alpha*L)*... sin(beta*L2*alpham*W)16*alpham^2*alpha^3*beta^2*c3*c2*... exp(alpha*L)*cos(beta*L2*alpham*W)+8*alpham^4*c2*c4*exp(alpha*L)*... alpha*beta^2*sin(beta*L+2*alpham*W)+8*alpham^2*c2^2*alpha^4*... exp(alpha*L)^2*beta*sin(2*alpham*W)+8*c1*alpha^4*c4*beta^3*... exp(alpha*L)*cos(beta*L+2*alpham*W)+8*alpha*L*alpham^2*... sin(2*alpham*W)*beta^5*c3^2+64*alpha^3*L*alpham^2*sin(2*alpham*W)*nu*.. c2*c1*beta^3+8*alpham^5*c2^2*exp(alpha*L)^2*beta*alpha^2*W4*alpham^2*... c3*c4*beta^2*alpha^3*sin(2*beta*L2*alpham*W)+beta^6*c4^2*alpha*... cos(2*beta*L2*alpham*W)8*c1*alpha^4*c4*beta^3*exp(alpha*L)*... cos(beta*L2*alpham*W)+32*alpham^3*alpha*beta^4*nu*c1*c3*... exp(alpha*L)*sin(beta*L)*W+32*alpham^5*c2*c4*exp(alpha*L)*alpha^2*... beta*sin(beta*L)*W+32*alpham^3*alpha^2*beta^3*nu*c1*c3*... exp(alpha*L)*cos(beta*L)*W+16*alpham^3*alpha^3*beta^2*nu*c3^2*... sin(2*beta*L)*W32*alpham^5*c2*c4*exp(alpha*L)*alpha*beta^2*... cos(beta*L)*W8*alpham^2*alpha^2*beta^3*c1^2*exp(alpha*L)^2*... sin(2*alpham*W)+16*sin(2*alpham*W)*c2*alpha^3*c4*beta^432*alpham^3*... alpha^4*beta*nu*c2^2*exp(alpha*L)^2*W+32*alpham^3*alpha^3*beta^2*nu*... c3*c2*exp(alpha*L)*sin(beta*L)*W64*alpham^3*alpha^2*beta^3*c1*c3*... exp(alpha*L)*cos(beta*L)*W8*alpham^2*alpha*beta^4*nu*c2*c3*... exp(alpha*L)*cos(beta*L2*alpham*W)+c4^2*alpham^4*alpha*beta^2*... cos(2*beta*L2*alpham*W)+8*alpham^2*alpha*beta^4*nu*c2*c3*... exp(alpha*L)*cos(beta*L+2*alpham*W)8*alpham^2*c3*c4*beta^2*alpha^3*... sin(2*alpham*W)4*c1^2*alpha^6*sin(2*alpham*W)*beta+32*alpham^5*c2*... c3*exp(alpha*L)*alpha*beta^2*sin(beta*L)*W8*alpham^2*c3*c4*beta^4*... alpha*sin(2*alpham*W)+8*alpham^2*c2^2*alpha^2*exp(alpha*L)^2*beta^3*... sin(2*alpham*W)+8*alpham^2*alpha^4*beta*nu*c3*c2*exp(alpha*L)*... sin(beta*L+2*alpham*W)+8*alpham^4*c1*c4*exp(alpha*L)*alpha*beta^2*... sin(beta*L+2*alpham*W)+8*c2^2*alpha^6*exp(alpha*L)^2*beta*W*alpham+32*. .. alpha^3*L*alpham^5*c1*c2*W*beta4*alpham^2*c3*c4*beta^4*alpha*... sin(2*beta*L2*alpham*W)8*alpham^4*c1*c3*exp(alpha*L)*alpha*... beta^2*cos(beta*L2*alpham*W)8*alpham^2*alpha*beta^4*nu*c1*c4*... exp(alpha*L)*sin(beta*L2*alpham*W)8*alpham^5*c1^2*exp(alpha*L)^2*... beta^3*W16*alpham^3*alpha^4*beta*c1^2*exp(alpha*L)^2*W32*alpham^3*. .. alpha^2*beta^3*nu*c2*c4*exp(alpha*L)*sin(beta*L)*W8*c2*alpha^3*c3*... PAGE 97 90 Appendix G Continued beta^4*exp(alpha*L)*cos(beta*L+2*alpham*W)64*alpha^3*L*alpham^3*c1*... c2*W*beta^3+8*alpham^2*alpha*beta^4*nu*c2*c4*exp(alpha*L)*sin(beta*... L+2*alpham*W)+32*alpham^3*alpha*beta^4*nu*c2*c3*exp(alpha*L)*... sin(beta*L)*W+2*beta^4*c4^2*alpham^2*alpha*cos(2*beta*L2*alpham*W)+... 2*alpham^2*alpha*beta^4*c3^2*cos(2*beta*L+2*alpham*W)4*alpha*L*... alpham^4*sin(2*alpham*W)*c4^2*beta^3 16*alpha*L*alpham^4*... sin(2*alpham*W)*c2*c1*beta^3+8*alpha*L*alpham^2*sin(2*alpham*W)*... beta^5*c4^2+16*c1^2*alpha^4*alpham^3*W*beta16*alpham^2*alpha^3*... beta^2*c4*c2*exp(alpha*L)*sin(beta*L2*alpham*W)4*c4^2*alpham^5*... sin(2*beta*L)*alpha^3*W32*alpham^3*alpha^2*beta^3*nu*c2*c3*... exp(alpha*L)*cos(beta*L)*W8*alpham^4*c1*c4*exp(alpha*L)*alpha^2*... beta*cos(beta*L+2*alpham*W)+8*c1*alpha^3*c4*beta^4*exp(alpha*L)*... sin(beta*L2*alpham*W)16*alpha*L*alpham^2*sin(2*alpham*W)*c3^2*... beta^5*nu16*alpha*L*alpham^2*sin(2*alpham*W)*c4^2*beta^5*nu+8*... alpha^3*L*alpham^5*W*c3^2*beta8*alpham^4*c2*c3*exp(alpha*L)*alpha*... beta^2*cos(beta*L2*alpham*W)32*alpham^3*alpha*beta^4*nu*c2*c4*... exp(alpha*L)*cos(beta*L)*W32*alpham^3*alpha^4*beta*nu*c4*c2*... exp(alpha*L)*sin(beta*L)*W+16*alpham^2*alpha^2*beta^3*c2*c3*... exp(alpha*L)*sin(beta*L+2*alpham*W)8*alpham^3*alpha*beta^4*c3^2*... sin(2*beta*L)*W+16*alpham^2*alpha^2*beta^3*c1*c3*exp(alpha*L)*... sin(beta*L2*alpham*W)8*alpham^4*c2*c3*exp(alpha*L)*alpha^2*beta*... sin(beta*L+2*alpham*W)16*alpham^2*alpha^3*beta^2*c1*c4*... exp(alpha*L)*sin(beta*L2*alpham*W)+8*alpham^4*c1*c3*exp(alpha*L)*... alpha^2*beta*sin(beta*L+2*alpham*W)+24*alpham^2*alpha^3*beta^2*nu*c1*.. c4*exp(alpha*L)*sin(beta*L2*alpham*W)4*alpham^4*c3*c4*alpha^3*... sin(2*alpham*W)+2*c3*beta^4*c4*alpha^3*sin(2*beta*L+2*alpham*W)4*... alpham^4*c1^2*sin(2*alpham*W)*beta*alpha^216*alpham^3*alpha^3*... beta^2*nu*c4^2*sin(2*beta*L)*W+16*alpham^3*alpha*beta^4*c3*c4*... cos(2*beta*L)*W+32*alpham^3*alpha^3*beta^2*nu*c3*c1*exp(alpha*L)*... sin(beta*L)*W+16*alpham^3*alpha*beta^4*nu*c3^2*sin(2*beta*L)*W+16*... alpham^2*alpha^2*beta^3*c2*c4*exp(alpha*L)*cos(beta*L2*alpham*W)+32*... alpham^3*alpha^2*beta^3*nu*c1*c4*exp(alpha*L)*sin(beta*L)*W32*... alpham^5*W*c2*c3*alpha^2*beta16*alpha^2*alpham^3*c2^2*W*beta^3c4^2*... alpham^4*alpha^3*cos(2*beta*L+2*alpham*W)c3^2*alpham^4*alpha^3*... cos(2*beta*L2*alpham*W)beta^6*c3^2*alpha*cos(2*beta*L2*alpham*W)... PAGE 98 91 Appendix G Continued beta^6*c4^2*alpha*cos(2*beta*L+2*alpham*W)+c4^2*alpham^4*alpha^3*... cos(2*beta*L2*alpham*W)+beta^6*c3^2*alpha*cos(2*beta*L+2*alpham*W)+... beta^4*c4^2*alpha^3*cos(2*beta*L2*alpham*W)+8*alpham^4*c1*c3*... exp(alpha*L)*alpha*beta^2*cos(beta*L+2*alpham*W)8*alpham^3*alpha^3*... beta^2*c3^2*sin(2*beta*L)*W+32*alpham^3*alpha^4*beta*nu*c1^2*... exp(alpha*L)^2*W4*alpham^4*c2^2*exp(alpha*L)^2*beta^3*... sin(2*alpham*W)+8*alpham^4*c1*c4*exp(alpha*L)*alpha^2*beta*... cos(beta*L2*alpham*W)32*alpham^3*alpha^3*beta^2*nu*c4*c2*... exp(alpha*L)*cos(beta*L)*W+8*alpha^3*L*W*c3^2*beta^5*alpham4*... alpham^4*c3*c4*alpha*beta^2*sin(2*alpham*W)2*alpham^4*c3*c4*alpha*... beta^2*sin(2*beta*L2*alpham*W)+2*alpham^4*c3*c4*alpha*beta^2*... sin(2*beta*L+2*alpham*W)4*alpha^3*L*alpham^4*sin(2*alpham*W)*c3^2*... beta4*L*sin(2*alpham*W)*beta^5*c4^2*alpha^34*L*sin(2*alpham*W)*... beta^5*c3^2*alpha^3+8*alpham^4*c2*c3*exp(alpha*L)*alpha^2*beta*... sin(beta*L2*alpham*W)+2*alpham^2*alpha^3*beta^2*c3^2*cos(2*beta*L+2*... alpham*W)+8*beta^4*c4^2*alpham^3*sin(2*beta*L)*alpha*W+8*alpham^2*... alpha^4*beta*nu*c4*c2*exp(alpha*L)*cos(beta*L2*alpham*W)+8*c1*... alpha^3*c3*beta^4*exp(alpha*L)*cos(beta*L2*alpham*W)+2*c3*beta^6*... c4*alpha*sin(2*beta*L+2*alpham*W)+4*c1^2*alpha^4*exp(alpha*L)^2*... beta^3*sin(2*alpham*W)24*alpham^2*alpha^2*beta^3*nu*c1*c4*... exp(alpha*L)*cos(beta*L+2*alpham*W)+16*alpham^3*c2^2*alpha^2*... exp(alpha*L)^2*beta^3*W16*alpham^2*alpha^3*beta^2*c3*c1*... exp(alpha*L)*cos(beta*L2*alpham*W)+32*alpham^3*alpha^2*beta^3*nu*... c1^2*exp(alpha*L)^2*W8*alpham^2*alpha^4*beta*c1^2*exp(alpha*L)^2*... sin(2*alpham*W)64*alpham^3*alpha^3*beta^2*c3*c1*exp(alpha*L)*... sin(beta*L)*W64*alpham^3*alpha^3*beta^2*c3*c2*exp(alpha*L)*... sin(beta*L)*W24*alpham^2*alpha^2*beta^3*nu*c2*c3*exp(alpha*L)*... sin(beta*L+2*alpham*W)+8*beta^2*c4^2*alpham^3*sin(2*beta*L)*alpha^3*... W beta^4*c4^2*alpha^3*cos(2*beta*L+2*alpham*W)16*alpham^4*... sin(2*alpham*W)*c1*alpha^2*c3*beta4*L*sin(2*alpham*W)*beta^7*c4^2*... alpha8*c1*alpha^4*c3*beta^3*exp(alpha*L)*sin(beta*L+2*alpham*W)... 16*alpha^3*L*alpham^4*sin(2*alpham*W)*c2*c1*beta4*alpha^3*L*... alpham^4*sin(2*alpham*W)*c4^2*beta+8*c2*alpha^4*c3*beta^3*... exp(alpha*L)*sin(beta*L+2*alpham*W)32*alpha^3*L*alpham^2*... sin(2*alpham*W)*c2*c1*beta^3+8*alpha^3*L*alpham^2*sin(2*alpham*W)*... beta^3*c4^216*alpha^3*L*alpham^2*sin(2*alpham*W)*c3^2*beta^3*... nu16*alpha^3*L*alpham^2*sin(2*alpham*W)*c4^2*beta^3*nu+32*c1*... alpha^4*c4*beta^3*exp(alpha*L)*sin(beta*L)*W*alpham8*c3*beta^6*c4*... cos(2*beta*L)*alpha*W*alpham4*c2^2*alpha^6*exp(alpha*L)^2*beta*... PAGE 99 92 Appendix G Continued sin(2*alpham*W)16*alpha^4*alpham^3*c2^2*W*beta16*alpham^2*... sin(2*alpham*W)*beta^4*c1*nu*c4*alpha+32*W*c2*alpha^4*c3*beta^3*... alpham32*c1*alpha^3*c3*beta^4*exp(alpha*L)*sin(beta*L)*W*alpham+8*... alpha^3*L*alpham^2*sin(2*alpham*W)*beta^3*c3^24*c3*beta^6*c4*alpha*... sin(2*alpham*W)4*beta^6*c4^2*sin(2*beta*L)*alpha*W*alpham32*c2*... alpha^4*c4*beta^3*exp(alpha*L)*sin(beta*L)*W*alpham+8*c2^2*alpha^4*... exp(alpha*L)^2*beta^3*W*alpham+32*c1*alpha^4*c3*beta^3*exp(alpha*L)*... cos(beta*L)*W*alpham+4*beta^4*c3^2*sin(2*beta*L)*alpha^3*W*alpham+8*... c2*alpha^4*c4*beta^3*exp(alpha*L)*cos(beta*L2*alpham*W)32*c2*... alpha^4*c3*beta^3*exp(alpha*L)*cos(beta*L)*W*alpham+4*beta^6*c3^2*... sin(2*beta*L)*alpha*W*alpham8*c1^2*alpha^4*exp(alpha*L)^2*beta^3*... W*alpham8*c3*beta^4*c4*cos(2*beta*L)*alpha^3*W*alpham+8*c1*alpha^4*... c3*beta^3*exp(alpha*L)*sin(beta*L2*alpham*W)8*c1*alpha^3*c4*... beta^4*exp(alpha*L)*sin(beta*L+2*alpham*W))/(alpha^2+beta^2)... /alpha/beta/alpham; function f=buckling(m,Nx,params) l=params.L; %Length of Gel w=params.W; %Width of gel alpham=m*pi/w; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); beta=sqrt(m^2*pi^2/(w^2)+sqrt(Nx/Def*m^2*pi^2/(w^2))); A=1; B=1; C=1; D=0; E=alpha; F=alpha; G=0; H=beta; I=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); J=(alpha^2*exp(alpha*l)alpham^2*nu*exp(alpha*l)); K=(beta^2*cos(beta*l)alpham^2*nu*cos(beta*l)); L=(beta^2*sin(beta*l)alpham^2*nu*sin(beta*l)); M=(alpha^3*exp(alpha*l)+alpham^2*(2nu)*alpha*exp(alpha*l)); N=(alpha^3*exp(alpha*l)alpham^2*(2nu)*alpha*exp(alpha*l)); O=(beta^3*sin(beta*l)+alpham^2*(2nu)*beta*sin(beta*l)); P=(beta^3*cos(beta*l)alpham^2*(2nu)*beta*cos(beta*l)); %Determinant of 4x4 Matrix f=A*(F*(K*PO*L)G*(J*PL*N)+H*(J*OK*N))... B*(E*(K*PO*L)G*(I*PL*M)+H*(I*OK*M))+... PAGE 100 93 Appendix G Continued C*(E*(J*PL*N)F*(I*PL*M)+H*(I*NJ*M))... D*(E*(J*OK*N)F*(I*OK*M)+G*(I*NJ*M)); function C=buckling2(m,Nx,params) L=params.L; %Length of Gel W=params.W; %Width of gel alpham=m*pi/W; nu=params.nu; %Poisson Ratio Def=params.D; %Flexural Rigidity alpha=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/Def*m^2*pi^2/(W^2))); beta=sqrt(m^2*pi^2/(W^2)+sqrt(Nx/Def*m^2*pi^2/(W^2))); A=[1 1 1 0; alpha alpha 0 beta;... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... (alpha^2*exp(alpha*L)alpham^2*nu*exp(alpha*L)) ... ( beta^2*cos(beta*L)alpham^2*nu*cos(beta*L)) ... ( beta^2*sin(beta*L)alpham^2*nu*sin(beta*L));... ( alpha^3*exp(alpha*L)+alpham^2*(2nu)*alpha*exp(alpha*L)) ... (alpha^3*exp(alpha*L)alpham^2*(2nu)*alpha*exp(alpha*L)) ... (beta^3*sin(beta*L)+alpham^2*(2nu)*beta*sin(beta*L)) ... ( beta^3*cos(beta*L)alpham^2*(2nu)*beta*cos(beta*L))]; c=[0;0;0;0]; Amat=[A c]; %Making C4=1 Aref=rref(Amat); C=[Aref(1,4);Aref(2,4);Aref(3,4);1]; 