MULTILAYER. B E A M ANALYSIS INCLUDING SHEAR AND G E O M E T R I C NONLINEAR E F F E C T S By Herman Ho Ming Kwan B.A.Sc, T H E UNIVERSITY OF BRITISH C O L U M B I A , 1985 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L E M N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E O F M A S T E R OF APPLIED SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES (CIVIL ENGINEERING D E P A R T M E N T ) We accept 1 his thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA May 1987 ©Herman Ho Ming Kwan, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia 1956 Main Mall Vancouver, Canada Department V6T 1Y3 DE-6C3/81) A bstract This thesis presents an analysis and experimental verification for a mul-tilayer beam in bending. The formulation of the theoretical analysis includes the combined effect of shear and geometric nonlinearity. From this formulation, a finite element program (CUBES) is developed. The experimental tests were done on multilayer, corrugated paper beams. Failure deflections and loads are thus obtained. The experimental results are reasonably predicted by the numerical results. Based upon this com-parison, a maximum compressive stress is determined for the tested beam. Finally, design curves for the tested beam are drawn using the deter-mined maximum compressive stress and the finite element program. ii C ontents page Abstract ii List of Tables v List of Figures vi Acknowledgement vii 1 Introduction 1 1 General Theoretical Analysis 6 2 Formulation of Theory 8 2.1 General Assumptions 8 2.2 Kinematic Relationships 9 2.3 Shear Strain Approximation 11 2.4 Stress-strain Relationship 13 2.5 Virtual Work Equation 14 2.6 Effect of Shear Stress Continuity Constraint 14 3 Finite Element Formulation 16 3.1 Interpolations 17 3.2 Virtual Work Equations 20 3.3 Newton-Raphson Method 23 3.4 Method of Computation 25 3.5 Numerical Integration 25 3.6 Consistent Load Vector 26 iii 4 P r o g r a m Ver i f i ca t ion - Compar i son w i t h P rev ious Resu l ts 28 4.1 Geometric Nonlinearity 28 4.2 Shear 32 II An Application: Multilayer Corrugated Paper Beam 35 5 Exper iment 37 5.1 Experimental Description and Results 37 5.2 Comparison of Experimental and Numerical Results 43 5.2.1 Data Consideration and Numerical Results 44 5.2.2 Discussion of Results 51 6 Des ign Curves 56 7 Conc lus ion 58 B ib l i og raphy 59 A A Tangent ia l Stiffness M a t r i x Te rm 61 B M a t e r i a l Non l inear i t y - Shear M o d u l i 63 B.l General Formulation 63 B.2 Piecewise Linear Shear Strain Assumption 63 B.3 Finite Element Formulation 66 B.4 Newton-Raphson Method 72 B.5 Determining the Unknowns 76 C Sample Ca l cu la t i on of Equivalent E last ic M o d u l i 78 D P r o g r a m L i s t i ng 81 iv List o f Tables 3.1 Locations and Weights of Integration Points 26 5.1 Thickness of Paper 39 5.2 List of Tests Performed 40 5.3 Tests Properties 40 5.4 Test Results of Loads and Midspan-Deflection at Failure 44 5.5 Experimental and Numerical Results 51 C. l Computed Equivalent Elastic Moduli 80 v List of Figures 1.1 Kao's Shear Strain Approximation 2 1.2 Cross-sectional Displacement for Three Bending Theories 3 2.1 Side View of Typical Beam Layout 9 2.2 Deflection of Beam 10 3.1 n-th Finite Element in x-coordinate 16 4.1 Comparison of Simply Supported Beam Problem 29 4.2 Comparison of Fixed Ends Beam Problem 30 4.3 Comparison of Buckling Problem 31 4.4 Comparison of Cantilever-Shear Problem 32 4.5 Comparison of Three-Span Beam Problem 33 5.1 Experimental Setup 38 5.2 Direction of Corrugation 38 5.3 Front View of Beam Tested in Machine Direction 41 5.4 Cross-Section View of Beam Tested in Machine Direction 42 5.5 Typical Load-Deflection Curve of Beam Tested in Machine Direction 43 5.6 Compression Crease Failure 44 5.7 Sinusoidal Shape of Corrugation 45 5.8 Geometric Properties of the Approximated Triangular Shape . . . . 46 5.9 Comparison of the Approximate Shapes of Corrugation 47 5.10 Comparison of Experimental and Numerical Results in Machine Di-rection 52 5.11 Comparison of Computed and Experimental Load-Deflection Curve in Cross-Machine Direction 54 5.12 Typical Shear Behaviour of Testing Direction 55 6.1 Loads Interaction Curves at Failure Stress of 8.96 MPa 57 C. l Distances d's of the 18-layers beam 79 vi A cknow ledgement I would like to express my gratitude to my advisor, Dr. R. 0 . Foschi, for providing me with his invaluable assistance and guidance throughout the duration of the research period and the writing of this thesis. The financial support from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Also, the testing materials supplied by MacMil lan Bloedel Limited is received with thanks. Finally, I would like to thank the U.B.C. Civi l Engineering Department staff, my friends and my parents for their assistance, advice, and encouragement in com-pleting this thesis. vii Chapter 1 Intro duction Beam bending analysis often neglects shear strains by assuming negligible shear contribution to deflection. However, this assumption is not always correct. For example, shear effects can be very important in bending of short, stocky members, glulam beams in high humidity environment, and multilayer beams with relatively weak layers in shear. In addition, the geometric nonlinear effect (large deformation) in beam bending is important for the study of beam behaviour under the interaction of axial and lat-eral loads. Elastic deformation of the beam under such loading conditions requires that large displacement terms be added to the usual small deformation equations. Beam analysis considering the effects of shear and large deformations is required for many applications. Combined action of axial and lateral forces on laminated beams, multilayer composite beams, or even corrugated cardboard beams are ex-amples to which this analysis may be applied. Many researchers have investigated the bending of multilayer beams with core layers weak in shear. Kao and Ross (1968) considered the variation of total energy in obtaining a system of equations which describes the bending behaviour of a mul-tilayer beam. Shear strains from the core layers weak in shear, bending strains from the stiff face layers and in-plane displacements of the stiff layers were included in the strain energy calculation. Shear strain in each core is assumed to be constant 1 over the thickness of the core. However, no shear stress continuity was imposed at the interfaces of the layers. The shear strains of the cores due to the in-plane dis-placements, were approximated as the rate of change of the in-plane displacements at the midplanes of the two stiff layers in contact (7„ = ^j) as shown in Figure 1.1. Also, the normal stresses perpendicular to the span were ignored; thus, all points Nodel Node 2 A u -- tx, + ] - iii /'th stiff loyer yth core Al ~- {/vt)th stiff _ i layer Figure 1.1: Kao's Shear Strain Approximation (Khatua and Cheung, 1973) on the same cross-section were assumed to have the same lateral deflection. Khatua and Cheung (1973) modified Kao's shear strain equation at each core to approximate the shear strains more accurately. The shear strains of the cores due to the in-plane displacements were then approximated as the rate of change of the in-plane displacements across the thickness of the weak cores. This different assumption was represented by a factor multiplied to Kao's shear strain equations. This factor was dependent on the thickness of the core and the adjacent stiff lay-ers. In addition, a finite element approach was taken to formulate an approximate solution to the problem. Further investigation by Foschi (1973) suggested that shear deformation in the stiff layers should be included in the analysis. Also, shear stress continuity should be imposed at the interfaces of the layers. Foschi proposed a piecewise shear strain 2 function to model the shear strains over the thickness of the entire plate. Shear strains in the stiff layers were represented by a quadratic polynomial over the thick-ness of the layers; while, constant shear strains were assumed in the weak cores. However, in-plane displacements of the neutral plane (2 = 0 plane) were neglected. Therefore, in the finite element formulation, the unknowns at each node were only the lateral deflection w, the slope tu', and the shear strain 7 at the midplane of each layer. In addition, symmetric cross-sections and small deformation were assumed in the analysis. Putcha and Reddy (1986) proposed a higher order plate theory to account for the additional deformation contributed by shear strains. Thus, the normals to the midplane (z = 0 plane) before deformation were now neither straight nor normal to the midplanes. This higher order theory is advantageous over the Reissner-Mindlin first-order theory (Figure 1.2) because this theory satisfied the zero shear stresses conditions required on the top and the bottom faces of the plate. Also, the shear f >x (c) HIGHER ORDER THEORY UNOEFORMEO (a) CLASSICAL KIRCHH0FF THEORY (b) MINOLIN-REISSNER THEORY Figure 1.2: Cross-sectional Displacement for Three Bending Theories (Ren, 1986) 3 correction coefficient required in Reissner's theory was no longer necessary. In addition, geometric nonlinear effect was included in the strain equations. Finally, a mixed finite element model was formulated to approximate the theoretical solution. Each node of the element contained eleven degrees-of-freedom: 1. the usual u, v and tw displacements, 2. rotations of the normals in the xz and yz planes (6x,0y) and 3. six moment resultants. Ren and Hinton (1986) later modified Reddy's higher order theory to develop a new finite element to investigate the simpler problem of laminated plate without geometric nonlinearity. Ren replaced the the rotation of the normals' degrees-of-freedom with the shear strains at the z = 0 plane. The u and v displacement functions were modified to contain the unknown shear strains. However, the zero shear stress conditions were still maintained. Di Sciuva (1986) proposed a'zig-zag displacement model' to account for shear effect on plate bending. This model incorporated a piecewise linear distribution of the in-plane, u and v, displacements. Shear stress continuity was imposed at the layers' interfaces. However, since the in-plane displacements were only linear, the shear strains in each layers were forced to be constant across the thickness of the layers. Thus, the model could not satisfy the zero shear stresses conditions on the top and the bottom faces of the plate. Therefore, shear stresses could not be directly obtained from the displacement model; instead, membrane stresses obtained from the model were substituted into the elasticity equilibrium equations to determine the shear stresses. The analysis also included the geometric nonlinearity terms in the strain equations. Finally, the model assumed that the plate's cross-sections were symmetric about the midplanes. This thesis utilizes many of the ideas mentioned in the above review to study the behaviour of multilayer beams in bending. These ideas include: 1. Strain energy contribution from bending and shear from all layers, 4 2. Finite element method to formulate the approximate solution to the problem, 3. Negligible normal stresses perpendicular to the span, 4. Shear stress continuity between each layer, 5. Inclusion of geometric nonlinear terms for combined lateral and axial loads, 6. Piecewise linear shear strains across beam depth, 7. Zero shear stresses conditions on the top and bottom faces of the beam. The finite element program CUBES is applied to the case of multilayer beams manufactured from corrugated paper. The theoretical solutions are compared to experimental results from bending tests of such beams. The program is also applied to the developement of strength interactions criterion when these corrugated beams are subjected to combined lateral and axial loads. Finally, a formulation to include material nonlinearity in shear is proposed. 5 Part I General Theoretical Analysis 6 This section discusses the general formulation of the analysis. Chapters 2, 3, and 4 are included in this section. Chapter 2 presents the general formulation of the beam bending theory which includes shear effect, geometric nonlinear effect, and multilayer beams. Chapter 3 describes a finite element formulation to approximate a solution to the theoretical problem. Finally, chapter 4 presents several compar-isons of results from the finite element program CUBES and other theoretical and numerical analyses. 7 Chapter 2 Formulation of Theory A virtual work approach is taken to analyze the problem of bending of a mul-tilayer beam. The general assumptions made in the analysis are first outlined. Kinematic relationships for strains and displacements are then developed. Finally, the governing equations are derived by applying the principle of virtual work. 2.1 General Assumptions Several assumptions are made to simplify the analysis: 1. Small strains are assumed. 2. Normal stresses perpendicular to the beam span are ignored; hence, across the beam depth, the lateral deflection w is assumed to be the same for all layers. 3. Elastic material properties are assumed. 4. Homogeneous material properties are assumed for each layer. 5. Solid rectangular sections are assumed for all layers. 6. A l l layers in the beam are in constant contact with each other, thus no dis-continuity exists between layers. 8 7. Out-of-plane warping of the beam is prevented. 8. Poisson effects are ignored. A typical beam cross-section with the layers' number is shown in Figure 2.1. 1 Gn-| G„ x — J Tn.." u = -\' . NL Figure 2.1: Side View of Typical Beam Layout 2.2 Kinematic Relationships Beam bending analysis begins with the consideration of strain. The strain equations can be obtained by looking at a point P and its deflection in the beam (Figure 2.2). Let us define a coordinate system [x,y,z) with origin at point 0 . The x-axis is parallel to the span and the plane y - z is the plane of the beam cross-section at O. Consider now a point P on this cross-section. The vertical distance between points O and P is then equal to z. Also, the displacements u and w are defined to 9 Figure 2.2: Deflection of Beam correspond with the directions x and z. The lateral deflections w(x) are assumed to be the same for all points on the same y — z plane. The axial displacement on the neutral axis is defined as u(x). Referring to Figure 2.2, if shear deformation is neglected, point P will deflect to P'. However, if shear is included, a term u*(x,z) must be added to the small deformation theory's axial displacement equation. The final axial displacement u(x,z) is thus expressed as u(x, z) = u(i) — z——h «*(x, z) (2-1) dx The strains at any point are defined by the equations _ du ^ 1 jdwy, ^ dx 2 dx ^ _dw ^du^^du^ dx dz dx 10 The last term in equation 2.2 is included because of the geometric nonlinearity (large deformation) consideration for the displacement w. The last term in equation 2.3 is approximated as ^ because j£ is approximated to be much smaller than 1. Substituting equation 2.1 into equations 2.2 and 2.3, the final strain equations become _ du drw ^ du* ^ 1 ^ du)y ^ 4) dx dx* dx 2 dx dw dw du* du* . . dx dx dz dz 2.3 Shear Strain Approximation Given the layout of a NL-layers beam shown in Figure 2.1, it will be assumed that the shear strain 7 I Z varies linearly over the thickness of each layer. Thus, two shear strain parameters per layer are needed to fully describe the shear strain distribution. The shear strain equation for the n-th layer can thus be expressed as ^ ) = ^ ^ + 7 ^ (2.6) It should be noted that 7 n and 7„ 0 are functions of the x-coordinate. The parameter c is a non-dimensional local coordinate in the z direction. At the point of maximum z, f equals to 1 in the layer while at minimum z, f equals to -1. In order to approximate the behaviour of continuous shear stress across the layers, the assumption of shear stress continuity is imposed. Thus, in addition to the virtual work requirement, the solution to the bending problem must satisfy the shear stress continuity constraint (see section2.6). Now, with given shear moduli Gn and Gn-i for n-th and (n-l)-th layers, the assumption of shear stress continuity between layers gives < ? n - l 7 n - l = G n 7 n 0 (2.7) Substituing 7„ 0 from equation 2.7 into equation 2.6, the shear strain in the n-th 11 layer is finally expressed as 7»U) - 7 n " + -Q—ln-l ^ (2.8J Integrating the above equation with respect to z gives the shear distortion equation u*(x,z) at any point (x,z) in the beam. Thus, u*(x,z) is given by u*(x,z) — f ^(x,s)ds = j i(x,s)ds— [ i(x,s)ds (2.9) JO J-h/2 J-h/2 Now, the global z coordinate in the n-th layer is transformed to the local coordinate (. For the n-th layer, zn is defined as the z value of the midplane (f = 0) and tn is defined as the thickness. The following transformation equation is obtained. z = zn + ^ (2.10) Equation 2.10 is then substituted into equation 2.9 which results in the following expression for the n-th layer's shear distortion. «»(*> ?) = £ [ ! f \ 7,(x, <r)rf<r] + y in{x, e)d<r K . . rLOCAL The number L O C A L is defined as the local coordinate value in the NA-th layer where the plane z = 0 is located; thus, NA is defined as the layer number for the point O. The last two terms in equation 2.11 are constants in z and are determined by the location of the z= 0 plane. But these two terms are functions of the x-coordinate. Finally, using equation 2.8, equation 2.11 can be changed to give the 12 shear distortion u*(x, z) at any point (x, z) in the n-th layer of the beam as n NA < ( * > * ) = £ 7 < ( x ) F ( n , 0 - £ 7 . ( * ) * , > , 0 (2.12) where = I + & - - ••»+i^ r'i - T +1 ' A <" -for » = 1,2,..., n — 1 F(n,n) = ^(i+Sl + h for t = n Ji"(JVX ,0 = | + M ± i ( l - A(ATA - ,•)) for t = 1,2,...,NA - 1 for t = JVA and A(/) = 0 if / > 1 and A(/) = 1 if / = 1 where / is the argument of the Afunction 2.4 Stress-strain Relationship The material properties expressed in terms of E and G are assumed to be elastic for each layer. Poisson effects are ignored. The bending and shear stress-strain relationships in the n-th layer can be expressed as °xn - EXneXn 13 (2.13) (2.14) or EXn 0 K M [ 0 GXZn J \ l x Z n f = [D(n)]{e(n)} (2.15) Another nonlinear effect not considered in the formulation of this analysis is the material nonlinearity. However, this effect can be important if the stress-strain curve deviates greatly from the linear approximation during deformation. This topic, specifically for the shear moduli, is discussed in Appendix B. 2.5 V i r t u a l Work Equation Knowing the stress-strain relationship, we must now determine a governing equation in order to find the unknown deflections. Applying a virtual displacement {6a} to the system, the resulting external and internal work (SW and SU ) done by the forces and the stresses in the system are respectively given as SU = jv{6e}T{a}dV SW = Jv{6a}T{F}dV (2.16) Equating the external and internal work done by the system gives the virtual work equation of Jv{6e}T{a}dV = Jv{6a}T{F}dV (2.17) 2.6 Effect of Shear Stress Continuity Constraint In addition to the virtual work requirement of equating the external and internal work done, a constraint of shear stress continuity between layers was imposed. Thus, 14 a 'constrained' virtual work equation would be derived which would satisfy the requirements of equating the work done and imposing the shear stress continuity. Such a constrained set of equations normally complicate the solution as the number of unknown parameters usually increases (Zienkiewicz, 1979). However, this shear stress continuity constraint actually simplifies the solution by reducing the number of unknowns in the problem. Because the set of shear stress continuity equations are simple linear equations, the constraint can be directly substituted into the virtual work equation. The substitution of the shear stress equations would thus eliminate the shear stress continuity constraint. Assumed displacement interpolations can now be applied to solve the bending problem. In our case, the unconstrained virtual work equation would contain two shear strain parameters from each layer in addition to the displacement parameters of u,w, and their derivatives. By applying the shear stress constraint as shown in equation 2.7, the displacement parameters of shear strain from each layer would be reduced from two to one parameter per layer in the 'constrained' virtual work equa-tion. This equation can now be solved using finite element method with assumed interpolations for the unknown displacements. 15 Chapter 3 Finite Element Formulation The finite element method is used to obtain approximate, numerical solutions to the set of governing equations derived from the concept of virtual work. A beam element with two end nodes is used in the formulation. Local coordinate £ is used in each element (Figure 3.1) along the x-axis. Each element has two end nodes, n and n+1, separated by the length Ax. The x-coordinate at £ equals to -1 of the n-th element is defined as x n . Ax Figure 3.1: n-th Finite Element in x-coordinate Thus, the x-coordinate of any point in the n-th element can be expressed as node n node rt + \ Element n x = x n + — (£ + 1) (3.1) and the differential, dx is (3.2) therefore dx Ax 2 (3.3) 16 The displacements u and tw and their 1-st derivatives u'(= j£) and w'(= *jr) are specified as the nodal degrees of freedom (DOF) for the interpolations. In addition to the four DOFs mentioned above, the nodal displacement vector for each node contains every layer's shear strain parameter at the node. Therefore, for a NL-layers beam, each node will have (4+NL) degree of freedoms per node. The degree of freedom of each element can then be assembled into a column vector called the displacement vector a. This vector takes the form: < l2n lNLn K+l un+l K+l l lNLn+l (3.4) 3.1 Interpolations Complete cubic interpolations are used to approximate the displacements u and tw within any element. It should be noted that, according to the compatibility requirement, displacement u needs only be a linear interpolation. However, during program implementation, a cubic interpolation was found to give a much improved approximation of the axial stresses. A linear interpolation is used to approximate the shear strain along x. Also, as previously described, the shear strain in each layer is approximated with a linear interpolation along z. For a complete cubic interpolation, four constant parameters are needed to define 17 the function. The displacement and its 1-st derivative at the two nodes provide sufficient parameteres to fully describe a cubic interpolation. The displacements un and wn can thus be expressed as Linear shear interpolations in both x and z direction are similar in form. For the x direction, the interpolation is 7,„(£,f) = : ^ ( l - 0 + ^ ^ ^ ( 1 + 0 (3.7) where 7,-„ is shear strain at the i-th layer of the n-th node and 7,-n+1 is shear strain at the i-th layer of the (n+l)-th node. For the z direction, the interpolation has already taken form in the previous section of shear strain approximation (eq. 2.6). The generalized displacements (u, w, and 7) can now be represented in vectorial form. The u-displacement can be written as a function of {a} in the following vectorial form «(0 = {^ l ( 0 } T {a} (3.8) where the vector {JVl(f)} is 18 {Nl(t)} = 0 A x l - 3 ( ^ ) ! + 2(i±i)'] 0 • , ( t t l ) ' _ a ( t ± l ) ' A x 0 Similarily, the displacement w can be expressed as term 1 st 2 nd 3 rd 4 t/i 5 th NL + 6th NL + Tth NL + 8th NL + 9th 2NL + 8th (3.9) (3.10) where {M(£)} is written as {M(£)} = A x [ l - 3 ( ^ ) 2 + 2 ( ^ ) ! (£±i)-2f«±i)1+r£tiv v 2 ; o o V 2 y 0 A x 0 0 term 1 st 2 nd 3 rd 4 r/i NL + 4 th NL + 5th NL + 6th NL + 7th NL + 8th 2NL + 8th (3.11) Finally, the linear shear interpolation in £ can also be expressed as a product of 19 vectors with the vector {iV3,(£)} equals to (3.12) 0 0 0 ( ¥ ) 0 0 0 term 1 st 2 nd i + 3th i + 4th i + 5th i + 6th NL + i + 7 th NL + i + 8th (3.13) 0 > 2NL + 8th In addition, derivatives of the displacements can be similarily expressed in the vector form. However, noting that the strain equations require derivatives with respect to the global coordinates x and z, the following equations are necessary to relate the derivatives in different coordinates. du* dz du dx dw dx du* 2 du 2 a*£ Ax dw 2 d£ Ax 3.2 Virtual Work Equations (3.14) Given equations 3.4 to 3.14, the strain equations (eqs. 2.4 and 2.5) for the n-th layer can now be expressed as 20 2 d{Nl} A x d£ - ( - S O 4 d 2 {M} A x 2 <f£2 1=1 1 4 + -2 A x • d{M} d{M} di t=i r {a} {KX{t ()f {a} + L-{a)T [MX(Z, f)\ {a} (3.15) = { m(e,f )} r {a} (3.16) Therefore, the strain equations are now expressed as functions of the displacement vector {a} and each strain contains two components: linear and nonlinear terms. Also, the strain vector {e} can be defined as Ui) \ where [5o(n)] is the linear component [flo(n)] = and [fii(n)] is the nonlinear component [B0(n)]{a} + [fl1(n)]{a} {KX}T {KXZ}T \{a}T\MX\ (3.17) (3.18) (3.19) Referring back to equation 2.17 of the virtual work equation section, (section 2.5) it can be seen that virtual strain equations are now needed in setting up the system of equations to be solved. Applying a virtual displacement to the strain equations, the virtual strain equations becomes we. ?)} = {lZti% }= [Bo{n)] {6a} + {B2{n)] {6a} (3-20) 21 where (S0(n)] is the previously defined linear component and [2?2(n)] is a nonlinear term which is written as = 2[Bi(n)] (3.21) Combining equations 2.15, 2.17, 3.17, and 3.20, the virtual energy equation can now be expressed as Sn = {6a}T fv([[B0(n)f + [B2(n)f] [D(n)} [[B0(n)} + [*,(»)] ] {a} - {F}) dV (3.22) Taking out the virtual displacement {6a} and setting the above equation to zero ( 611 = 0 ) , the equation takes the final form of {F} = [fv[[Bo(n)]T + [B2(n)}T][D(n)} [[B0(n)) + [B1(n)]]dV] {a} [K(n)}a (3.23) where the right hand side of the equation excluding {a} can be symbolized as [*(»)] = I [[BoHf + [52(n)]T] [D{n)\ [[BQ(n)} + [B^n)} }dV (3.24) with [-fiT(n)] being defined as the elemental structural stiffness matrix for the n-th layer. Expressing in local coordinate £ and (, the same matrix becomes \K(n)} = | ^ A j / £ df £ [[B0(n)]T + [B2(n)]T] [/>(»)] [[B0[n)) +[B1{n)]]dt (3.25) where Ay is the beam width. This matrix is determined for every layer in each element. The elemental structural stiffness matrices are then assembled into the 22 [B2(n)\ = {a}T[MX} {oy global structural stiffness matrix by summing across the beam depth and span. Because second order effect of geometric nonlinearity is included in the analysis, an iteration scheme is needed to solve the system of equations. The standard Newton-Raphson method discussed in the next section is used to solve the problem. 3.3 Newton-Raphson Method The Newton-Raphson method is a commonly used technique to solve nonlinear equations (Zienkiewicz, 1979). This method uses a first order approximation tech-nique to solve the equations through iterations. The first order approximation is stated as d$(a)' d{a} {Aa} (3.26) {$(o + Aa)} = {$(a)} + where {$} is a function of the displacement vector {a} . Equation 3.26 can be re-arranged to become W*(a)} 1 _ 1 d{a} (3.27) {Aa} = {{$(a + A a ) } - { $ ( a ) } } The value {Aa} can then be compared against an acceptable tolerance to determine whether further iteration is needed to obtain a sufficiently accurate solution {a}. From equation 3.27, it is obvious that [ d j * ^ ] n a s t° D e determined before the Newton-Raphson method is used. In our case, let {611} = {$(a)} = [fy [[B0(n)]T + [B2(n))T] [D(n)) x [ [B0(n)} + [B x(n)] }dv] {a} - {F} (3.28) where {$(a)} is now defined as the column vector {511}. The solution to the finite element problem is finding the displacement vector {a} which results in {$(a)} equal to or near zero. However, since the equations are nonlinear, a direct solution to the equations was not possible. 23 Remembering that the matrix [Bo] and forces {F} are not functions of {a}, the 1-st derivative can be obtained by differentiating equation 3.28 with respect to {a} and using equationd 3.20. The resulting equation is d{$(a)} d{a} Also, from the stress equation (eq. 2.15), the above term of is determined as d{a} d{a} (3.30) In addition, after some special manipulation (see Appendix A), the term dj^| T {a} is found to equal 'd[BjT d{a} W Ax MMX] (3.31) for our beam element. Substituting equations 3.30 and 3.31 into equation 3.29, d^jSaft then takes the final form of d{a} = j T ax\MX] + [[B0)T + [B2)T] [D][[B0] + [B2]]dV = [Kt] (3.32) Finally, this term [^j^] is called the tangential stiffness matrix [Kt]. Noting that both terms on the right-hand-side of equation 3.32 are symmetric matrices, computer storage can be minimized by storing the matrices in vector form. With the matrix [Kt] known, the following solution scheme is used to determine the approximate solution vector {am+i} in the (m+l)-th iteration. {*m(a)} + [Ktm]({am+1}-{am}) = {0} (3.33) {Ktm]{am+i} = [Ktm]{am}-{*m(a)} {am+1} = {am} - [KtS1 {*m(a)} The solution scheme is repeated until the displacement vector {Aa} < {tol}. 24 3.4 Method of Computation The Cholesky decomposition method of solution is used to determine the vector {Aa}. Symmetry and bandwidth are considered in storing the matrices. Before applying the Cholesky method, boundary conditions are first applied to [Kt] and { $ ( 0 ) } . For zero displacement boundary conditions, zeros are placed into the off-diagonal row and column of the specified degree of freedom in [Kt] while a zero is placed into the same DOF of the {$(a)} vector. Also, the value 1 is placed into the diagonal term of the specified DOF in the [Kt] matrix. Finally, the following procedures from the Cholesky method are used to solve the system of equations. 1. Decompose [Kt] 2. Solve {Aa} = [Kt]~l {$(a)> 3.5 Numerical Integration Since the integrals in the previous sections are very complicated, closed form solution is difficult to obtain; therefore, numerical integration is used. Gaussian quadrature scheme is applied because the scheme is most suitable with the local coordinate variation of -1 to 1. Full integration is used to compute the integrals. The maximum order of polynomial appearing in the integrals will determine the number of Gaussian points necessary to accurately integrate the function. The term [B2] contains a 4-th order polynomial in £ and [B2] is then squared in the [Kt] matrix. Thus, the highest order polynomial in the integrals is of 8-th order. Knowing that a n points Gaussian integration will integrate exactly a (2n-l)-th order polynomial, a 5-point Gaussian scheme is used for £ in the numerical integration. Following the same procedure, a 3-points Gaussian scheme is determined to be necessary for the f coordinate. Therefore, the integral can be represented as 25 point locations weights 5-points 1 integration 2 3 4 5 -0.9061798459 -0.5384693101 0.0000000000 0.5384693101 0.9061798459 0.2369268850 0.4786286705 0.5688888889 0.4786286705 0.2369268850 3-points 1 2 3 -0.7745966692 0.0000000000 0.7745966692 0.5555555556 0.8888888889 0.5555555556 Table 3.1: Locations and Weights of Integration Points /l -1 n m / F(e,fm = EE^ ( 6.&W (3-34) • 1 J - 1 t=ij=i where F(£,f) is any function in the coordinates (£, f). The locations and weights for a 5-point and a 3-point Gaussian scheme are given in Table 3.1. Since numerical integration is used, bending and shear stresses are computed only at the Gaussian points of each element. If stresses at any other point are desired, approximation such as stresses interpolation can be used. 3.6 Consistent Load Vector Because finite element method is used in the analysis, consistent load vector must be used if the applied load is to be represented exactly. The consistent load representation for a concentrated load is simply the load value placed in the specific degree of freedom. However, the consistent load representation for a distributed load is somewhat more complicated. The consistent load vector is represented by the following equation {Q} = ( ? ( 0 W 6 K (3-35) where g(f) is the load function and {M( f )} is the shape function within an element. In this program, only uniformly and linearly distributed loads within an 26 element are accounted for. For example, assuming an element is subjected to a laterally distributed load which varies linearly from qx to q2, the consistent load vector is determined by the following procedure. The load function g(£) within the element is (3.36) Substituting this equation into equation 3.35 and performing the integration, the consistent load vector becomes {Q} 20 " 1 20 A z 2 _ , Ax 2 „ "20"9l + -30-92 20 - A i 2 (3.37) 30 "9l 20 A simple check can be done on this equation by determining the better known, uniformly distributed case. Equating qi = q2, the consistent load vector of equa-tion 3.37 is {Q} = ( §Axgi 12 9l IT* (3.38) which corresponds with the uniformly distributed load case. 27 Chapter 4 Program Verification -Comparison with Previous Results The finite element program CUBES was developed based upon the formulation discussed in chapters 2 and 3. The program's results were compared to referenced results; the two sets of results were in close agreement with each other. For the geometric nonlinearity effect, two results given by Timoshenko and Woinowsky-Krieger (1959) and one given by Popov (1968) were compared to the program's numerical solutions. For the shear effect, Popov's theoretical solutions and Foschi's numerical solutions were compared to CUBES ' numerical results. 4.1 Geometric Nonlinearity Two beam problems from Timoshenko were used to verify the program's result. The first problem considered was a simply supported beam with supports at a fixed distance apart subjected to an uniformly distributed load. Maximum deflections of the beam under different load value were calculated using Timoshenko's theory and the program CUBES. The analytical and numerical results plotted in Figure 4.1 were almost identical. 28 UNIFORMLY LOADED SIMPLY-SUPPORTED BEAM WITH NO ROLLER E=206,850MN/m**2 L=1.27m. d=12.7mm. b=25.4mm. E E X a 30-i 25 20-O^r * i i < i ~r 500 ' 1000 1500 2000 k=qL**4/[Eld] — i — 2500 Legend TIMOSHENKO X CUBES program 1 3000 Figure 4.1: Comparison of Simply Supported Beam Problem The second problem of deflection of a fixed ends beam under uniformly dis-tributed load was also determined using Timoshenko's theory and the program CUBES. The results for this case are shown in Figure 4.2. The program's approxi-mate solutions were also in close agreement with Timoshenko's theoretical solutions. Finally, CUBES ' results for the buckling of a simply supported beam under combined lateral and axial loads were compared to Popov's theoretical solutions. The results are shown in Figure 4.3. 29 UNIFORMLY LOADED FIXED ENDS BEAM E=206,850 MN/m**2 l=1.27m. d=12.7mm. b=25.4mm. 2 5 -2 0 -£ 1 5 " £ X n q £ ^ 10-/ / / / 1 1 1 1 1 1 s / ' s / L 5 -n > Legend TIMOSHENKO X CUBES progrom u ^ c s 1 1 1 1 I ) 500 1000 1500 2000 2500 k=qL**4/[Eld] 3000 Figure 4.2: Comparison of Fixed Ends Beam Problem 30 E E X D E BUCKLING OF A SIMPLY SUPPORTED BEAM WITH ROLLER E=13,790 MN/m**2 l=1.524m. d=38.1mm. b=25.4mm. F=44.48 N Pcr=6859 N 1500 1000-500 -0.2 0.4 0.6 P/Pcr 0.8 —I 1.2 Legend CUBES linear u 10 ELEMENTS CUBES cubic u 10 ELEMENTS POPOV Figure 4.3: Comparison of Buck l ing P rob l em A s mentioned in chapter 3, a cubic u interpolation was chosen over a linear u interpolation in the program. For smal l number of elements used, the results of the linear u interpolation near cr i t ical load converge to the theoretical buckling solutions very slowly; thus, a cubic interpolation of u was used to improve the convergence. A s shown in Figure 4.3, this cubic u interpolation approximated the buckling results quite accurately even when small number of elements (ten elements along span) were used. A l l three of the above comparisons indicated that the program can accurately approximate the geometric nonlinear effect. 31 4.2 Shear For shear effect verification, a cantilever problem from Popov was used to compare with CUBES's results. Popov assumed a parabolic shear strain distribution across the beam depth. To include shear in CUBES's result, the beam was divided into equivalent fictitous layers thus approximating the shear strains with a piecewise linear distribution. The free end deflections were computed using two different approximating distribution: a 2-layers and a 4-layers approximation. It should be noted that, except for ux, all the degrees-of-freedom were specified as O's at the fixed end of the cantilever. Figure 4.4 shows the results of the two analyses. BENDING WITH SHEAR OF A CANTILEVER WITH CONCENTRATED LOAD P AT THE FREE END. E=206,850 MN/m**2 L=0.254m. d=25.4mm. b=25.4mm. 5 ELEMENTS 1.035 - i 32 The parameter a was defined as the ratio of total deflection with shear over flexural deflection. C U B E S ' predicted results were quite close to the theoretical results. Also, the piecewise shear strains distribution approaches the parabolic dis-tribution with greater number of fictitous layers. Figure 4.4 shows the comparison of the results. Finally, deflections, including shear contribution, of a three-span simply-supported beam were determined using the program CUBES. These results were compared and found to be near Foschi's results (1973) as shown in Figure 4.5. DEFLECTION' OF THREE-SPAN SIMPLY SUPPORTED BEAM WITH ROLLER, UNIFORMLY DISTRIBUTED LATERAL LOAD P COMPARISON OF FOSCHI'S RESULTS WITH CUBES - 4 ELEMENTS/SPAN Legend FOSCHI G2= 6.895 MN/m"2 FOSCHI C2= 68.95 m/m,,2 CUBES G2= •Foschi 6.895 MN/m**2 CUBES G2= 68.95 MN/m"2 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 SPAN L m. Figure 4.5: Comparison of Three-Span Beam Problem Four elements per span were used in the program to obtain C U B E S ' results shown in Figure 4.5. Only the linear results from the program C U B E S were used 33 in the comparison since Foschi's formulation did not include geometric nonlinear effect. The results from CUBES were not identical to Foschi's results because Foschi assumed a constant shear strain in the weak cores and quadratic shear strain in the stiff layers. This assumption was different with our present assumption of linear shear strains in all layers. Also, Foschi's formulation ignored u, the axial displacement of the z = 0 plane. The above two comparison with Popov and Foschi indicated that the program's shear approximations were also accurate. Therefore, combined with the geometric nonlinear approximations, the program could readily approximate a solution to the bending of a multilayer beam with shear and geometric nonlinear effects included. 34 Part II An Application: Multilayer Corrugated Paper Beam 35 This section focus on a specific application of the proposed theory: the bending of multilayer corrugated paper beams. Chapters 5 and 6 are included in this section. Chapter 5 describes the setup of a one-third point bending test of the paper beams and compares the experimental results to the program's numerical results. Chapter 6 shows a number of loads interaction design curves obtained from the program CUBES. 36 Chapter 5 Experiment A bending test was performed on multilayer corrugated beams. Results from this test were compared to the numerical results from the program CUBES. A strength failure criterion for the specific outer liner was then determined based upon the above comparison. 5.1 Experimental Description and Results A one-third point bending test was done to determine the midspan deflection of the tested beams. The beams were simply supported at both ends and were allowed to move axially (see Figure 5.1). Each beam was made from nine-layers of paper: one layer of 90-lb paper liner, four layers of 42-lb paper liner, and four layers of 26C paper corrugation. Each layer had orthotropic properties. The corrugation was machined in one direction only; thus, two directions were defined for the corrugated layer: 1. machine direction and 2. cross-machine direction. Machine direction (Figure 5.2) was the direction of the approximate sinusoidal waves while cross-machine direction was the direction perpendicular to the machine direction in the x — y plane. 37 The thickness of the paper was measured and listed in Table 5.1. The amplitude / and period Lc of the corrugation were respectively measured as 3.58mm. and 7.80mm.. Material Type Thickness (mm.) 90-lb 42-lb 26C 0.635 0.305 0.203 Table 5.1: Thickness of Paper The test was done in the U.B.C. materials lab on the SATEC machine. The relative humidity and temperature of the testing environment was measured daily and were respectively found to be around fifty percent (50%) and sixty-eight de-gree Fahrenheit (68°F). Also, before testing, each beam was placed overnight in the testing room to reach atmospheric equilibrium. The load was measured with a 444.8 N (100 lb.) load cell and applied at a displacement controlled rate of approx-imately 7.62mm./min.. This load cell accurately measured the applied loads of the tests which were in the 88.96-355.84 N (20-80 lb.) range. The midspan displacement was measured using a LVDT. Eighteen-layers beams were made from gluing two nine-layers boards together. After applying glue to both sides of the interface, the boards were pressed to-gether and left to bond overnight. Constant pressure of approximately 13.79 kN /m 2 (2.0 lb/in 2) was applied to the beam during the bonding period. The glued board was then checked to ensure that no crushing had occured from the applied pressure. Six different tests were done on the cardboard beams. They are listed in Ta-ble 5.2. The material properties of each of the tested beam are given in Table 5.3. 39 Test Test Layer Sample # Direction Arrangement Size 1 machine 90-26C-42 17 2 machine 90-26C-42-26C-90 11 3 machine 90-26C-42-90-26C-42 2 4 machine 90-26C-42 10 5 machine 42-26C-90 4 6 cross-machine 90-26C-42 2 Table 5.2: List of Tests Performed Test Beam No. of Beam Depth # Span (m.) Layers db (mm.) 1 0.6 9 16.28 2 0.6 18 32.56 3 0.6 18 32.56 4 0.3 9 16.28 5 0.3 9 16.28 6 0.9 9 16.28 Table 5.3: Tests Properties Figure 5.3 and 5.4 respectively show a front and a cross-section view of a 90-26C-42 beam tested in machine direction. The load-deflection curve was recorded during each test (Figure 5.5). Each beam was considered to have failed when a compressive crease had fully developed in the outer liner at the compression side of the beam (Figure 5.6). The complete devel-opment of this crease could also be seen in the load-deflection curve (see Figure 5.5) as a significant sudden decrease of the applied load . This definition of failure was applied to bending in both machine and cross-machine directions. The averages and standard deviations of the experimental failure loads and deflections for all the tests are given in Table 5.4. 40 41 Figure 5.4: Cross-Section View of Beam Tested in Machine Direction 42 f a i l u r e 0-18 Midspan D e f l e c t i o n Figure 5.5: Typical Load-Deflection Curve of Beam Tested in Machine Direction 5.2 Comparison of Experimental and Numerical Results The program's approximate solutions were compared to the experimental re-sults. However, certain data required special consideration because of the card-board's unique properties. The strength of the paper was highly dependent on the testing conditions of bending direction and relative humidity. Also, the program assumed a solid rectangular cross-section in each layer. This cross-section would provide a much higher bending stiffness than the actual cross-section in a corru-43 Figure 5.6: Compression Crease Failure Load Midspan Deflection Test # Average Standard Deviation Average Standard Deviation (N-) (N.) (mm.) (mm.) 1 69.4 8.4 12.45 1.40 2 162.6 8.2 6.10 0.58 3 158.4 0.7 6.86 0.25 4 68.3 7.1 4.83 0.51 5 45.6 1.8 3.81 0.13 6 40.8 0.0 64.77 1.98 Table 5.4: Test Results of Loads and Midspan-Deflection at Failure gated layer; thus, an equivalent elastic modulus was used to better approximate the beam's true bending stiffness. After the data were determined, results from the program were then compared to the experiemental results. 5.2.1 Data Consideration and Numerical Results Figure 5.7 shows the properties of the corrugation which has an assumed sinu-soidal shape. Special considerations were given to the following data: 1. Shear moduli of all the layers, 44 Figure 5.7: Sinusoidal Shape of Corrugation (Timoshenko,1959) 2. Elastic moduli of the liners, 3. Elastic moduli of the corrugated layers, The elastic and shear moduli values were given by MacMil lan Bloedel Research. The elastic moduli for each type of paper were determined from tension tests. The shear moduli for the corrugation were determined from direct shear tests of the complete corrugated sections. Shear Moduli The shear modulus Gxz was used in each layer for bending in the machine direction. Gyz was used for bending in the cross-machine direction. In addition, effect on the shear modulus due to relative humidity should be considered. Shear moduli at three different relative humidity values were given. The shear moduli corresponding to the measured relative humidity were then linearly interpolated from these three points. An alternate method of determining the shear modulus Gxz of the corrugation was done using the program NISA. Three approximate shapes of the corrugation were modelled in NISA to obtain shear stress-strain relationships in the x-z direc-tion. The three shapes are a straight-line, triangular shape, a sinusoidal shape, and a semi-circular shape. A single corrugation (see Figure 5.8) spanned over one wave 45 Figure 5.8: Geometric Properties of the Approximated Triangular Shape length Lc was used and the ends are assumed to be pin-ended with no roller. A concentrated force P was applied at the apex of the corrugation to produce the shearing action. The elastic modulus of the 26C paper was assumed to be constant. The shear strain was approximated as the displacement A over the amplitude /. The shear stress was approximated as the force per unit width P over the wave length Le. The resulting shear stress-strain relationships were plotted in Figure 5.9. The curve for the triangular shape is shown to be linear. But the curves for the sinusoidal and semi-circular shapes are shown to be nonlinear. The constant value of Gxz given by MacMil lan Bloedel Research falls within these two nonlinear curves. It is interesting to note that, although the elastic modulus of the paper was kept constant, the corrugation's geometry produced a nonlinear shear stress-strain relationship. Elastic M o d u l i of Liners For bending in the machine direction, elastic modulus Ex was used in each layer; while Ey was used for bending in the cross-machine direction. Relative humidity effect on the liners' elastic moduli was also considered using the same interpolation 46 COMPARISON OF CORRUGATION SHAPES Csl E 3 0 . 0 - , 2 5 . 0 2 0 . 0 8 UJ CC I— CO Q ; 1 0 . 0 < UJ co 5 . 0 0 . 0 0 . 0 0 0 . 0 2 r;„ ().(«):i»(ikN/iiiiii! 0 . 0 4 0 . 0 6 0 . 0 8 SHEAR STRAIN (unitless) o.io 0 . 1 2 Legend TRIANGULAR SINUSOIDAL SEMI-CIRCULAR GIVEN Gxz Figure 5.9: Comparison of the Approximate Shapes of Corrugation 47 method discussed in the shear moduli section. Elastic Moduli of Corrugated Layers In addition to the above considerations for bending direction and relative hu-midity, the elastic moduli of the corrugated layers should be modified because the program assumed a solid, rectangular cross-section in every layer. This assumption was not correct for corrugated layers. To correct for this assumption, equivalent elastic modulus for each corrugated layer should be computed and placed into the data file. To determine an equivalent elastic modulus, the layer's bending stiffness contribution in both x— and y—directions were equated for the corrugated and the solid layer: EXeqIXeq = ExIXcor(= Dx) (5.1) EyJVcq = EvIVcor(=Dy) (5.2) where EXeq, Eytq = equivalent elastic modulus for bending in the x— and y—directions, I x , Iycq = moment of inertia about the beam's global x— and y-axes of the assumed solid rectangular cross-sections, Ex,Ey = elastic moduli of the paper for bending in the x- and y-directions, IXeor, Iycor = moment of inertia about the beam's x- and y-axes of the corrugated layer. Dx,Dy = bending stiffness in the x-and y-directions If / = s(s is defined below), the corrugation would become a flat rectangular 48 cross-section of depth h. Its contribution to the beam's bending stiffness EI would be EI = E— + Ehd2 (5.3) 12 v ' On the other hand, if / = 0 the corrugation would not contribute to the beam's EI. It may be assumed (Timoshenko, 1959) that for the intermediate values of / (Figure 5.7), the bending stiffnesses in the x and y directions are: Dx, = - E x - (5.4) Dy, = EyI (5.5) where h is the thickness of the paper. The length of the arc for a half-wave is represented as s and is given as and 1 _ 0-81 l + 2.5(£)2 Because these bending stiffness were determined about the local axes x' and y' of the corrugation, the parallel axis theorem must be applied to determine the bending stiffness about the global axes x and y of the beam. Thus, if the distance in z direction between x' and x (also y' and y) was given as d, then the bending stiffness of the corrugation about the beam's axes (z and y) were obtained as Dx = Dx, + -Exhd2 , (5.6) Dy ' Dyl + -Eykd 49 (5.7) The second term in equation 5.6 was determined using the same linear approx-imation factor) that Timoshenko had presented (see equation 5.4). The second term in equation 5.7 was determined from finding the area per unit length of the cross-section (equal to y ) and then applying the parallel axis theorem. Applying the criterion in equation 5.2, the equivalent elastic moduli were then given as Applying equations 5.8 and 5.9 to the tested beams' corrugated layers, the equiv-alent elastic moduli Ex and EVeq were respectively found to be around 55.2 MPa and 262.0 M P a for all corrugated layers. The equivalent elastic moduli of all cor-rugations did not deviate greatly from the above two values (sample calculation in Appendix C). Numerical Results Numerical results of the tested beam at the experimental average failure loads were obtained from the program. Table 5.5 shows a comparison of the experimental average deflections and the computed deflections at the average failure loads. 50 Test # Sample Size Exper. Computed @ Failure Load Average Average Maximum Failure Failure Deflection Compressive Load, P (N) Defl. (mm.) (mm.) Stress (MPa) 1 17 69.26 12.45 13.63 10.1 2 11 162.35 6.10 6.71 8.0 3 2 158;13 6.86 6.82 8.0 4 10 68.23 4.83 4.46 9.1 5 4 45.59 3.81 2.98 10.2 6 2 40.03 64.77 40.08 6.4 Table 5.5: Experimental and Numerical Results *note: test # 6 was tested in the cross-machine direction 5.2.2 Discussion of Results From Table 5.5 above, the computed and experimental deflections in tests #1-4 were observed to be in close agreement with each other. Figure 5.10 shows a plot of the means and standard deviations of the experi-mental deflections vs. the computed failure deflections. The computed deflections of tests #1-4 are all located within the determined standard deviations. Also, Fig-, ure 5.5 shows that the load-deflection curves in tests #1-4 are almost linear. Based on the above comparison, the program may be said to provide a reasonable approximation to the test problem. Maximum compressive stresses in the 90-lb liner were then determined from the program for bending in the machine direction and the specific relative humidity in tests #1-4. By averaging the four tests' values, a maximum compressive stress of 8.96 MPa was obtained as the compressive failure stress of the 90-lb liner in the machine direction. This failure criterion was then imposed to determine load-interaction design curves for machine direction bending of any corrugated beam with the compressive crease failure occuring in the 90-lb liner. Result of test #5 was different because the compressive failure had occurred in the 42-lb liner. Thus, the maximum compressive failure stress was not included in 51 COMPARISON OF EXPERIMENTAL AND ANALYTICAL DEFLECTION RESULTS AT FAILURE LOAD Figure 5.10: Comparison of Experimental and Numerical Results in Machine Direction 52 determining the 8.96 MPa compressive stress value. Also, the computed deflection was significantly lower than the experimental average deflection. This difference may be attributed to nonlinear elastic moduli of the 90-lb liner in tension. In test #6, the beams were bent in the cross-machine direction. The experi-mental average deflection was much larger than the computed deflection. However, extrapolating a linear solution of 45.72 mm. from the experimental load-deflection curve (Figure 5.11), the computed deflection of 40.08 mm. would be a reasonable approximation to the linear solution. The concave shape of the experimental curve in Figure 5.11 indicated that the beam was softening during the test. Nonlinear effect had contributed significantly to the beam's deflection. This additional nonlinear deflection may be caused by nonlinear shear moduli of the corrugated layers. The deflections of the beams tested in the machine direction were in the small deformation range of ^ < 1 . Previous discussion of the shear moduli modelled by NISA has shown that the shear moduli in machine direction is highly nonlinear but the experimental curves are shown to be almost linear. Therefore, the beam's deflection must be occurring in the small shear strain range, thus producing a linear load-deflection curve. However, for beams tested in the cross-machine direction, the deflections were in the large deformation range of ^ >> 1. The above explanation of different deformation range can be observed in Figure 5.12. Typical shear stress-strain relationships for both test directions are assumed and the shear behaviour for the two directions are indicated. This difference may account for the linear and nonlinear experimental curves for the two different bending directions. 53 H l d s p a n D e f l e c t i o n ( I n . ) Figure 5.11: Comparison of Computed and Experimental Load-Deflection Curve in Cross-Machine Direction 54 SHEAR BEHAVIOUR OF TESTING DIRECTION SHEAR STRAIN Figure 5.12: Typical Shear Behaviour of Testing Direction 55 Chapter 6 Design Curves Design curves of load-interaction are produced by applying the maximum com-pressive failure stress criterion to the program's stress results. An example of load-interaction curve is presented for the bending of a simply supported beam subjected to an axial load and a linearly distributed lateral load The failure stress criterion of 8.96 MPa (1300 psi.) is applied to the outer liner (90-lb liner) on the compressive side of the beam. Combination of axial and lateral load values for spans of 1.219 m. (48 in.) and 1.829 m. (72 in.) are determined from the program. The results are plotted in Figure 6.1. Similar plots can be done to study the effects of different layer arrangements, elastic moduli, shear moduli, maximum failure stress value, etc., on the allowable load combination. 56 INTERACTION CURVE OF A SIMPLY SUPPORTED BEAM WITH ROLLER, 9 0 - 2 6 C - 4 2 - 4 2 - 2 6 C - 9 0 LAYER ARRANGEMENT MAX. ALLOW. COMP.=0.00896 kN/mm**2 W=25.4mm. 300-1 q, MAX. LATERAL LOAD VALUE (N/m.) Figure 6.1: Loads Interaction Curves at Failure Stress of 8.96 M P a 57 Chapter 7 C onclusion A finite element program has been developed using multilayer beam elements. Shear and large deformation effects are included in the analysis. The program is shown to accurately predict solutions which agree with several solutions found in the literature. Also, an experiment was done to measure the deflection of multilayer, corru-gated, cardboard beams which are weak in shear strength. The experimental results for bending in the machine direction were accurately predicted by the program's nu-merical results. A procedure to obtain combined axial and lateral loads interaction curves was presented using maximum compressive stress as a failure criterion. Further research relating to this subject should include: bending in the cross-machine direction, shear modulus nonlinearity, extension to plate bending, a shear stress approximation across the beam depth and dynamic loading. 58 B ibliography [1] DI SCIUVA, M. 1986. Bending, Vibration and Buckling of Simply Supported Thick Multilayered Orthotropic Plates: An Evaluation of a New Displacement Model. Journal of Sound and Vibrations 105(3), pp. 425-442. [2] FOSCHI, R.O. 1973. Deflection of Multilayer-Sandwich Beams with Applica-tion to Plywood Panels. Wood and Fiber, Vol. 5, No. 3, pp. 182-191. [3] K A O , J.S. and ROSS, R.J. 1968. Bending of Multilayer Sandwich Beams. A IAA Journal, Vol. 6, No. 8, pp. 1583-1585. [4] K H A T U A , T.P. and CHEUNG, Y.K. 1973. Bending and Vibration of Mult i-layer Sandwich Beams and Plates. International Journal for Numerical Meth-ods in Engineering, Vol. 6, pp. 11-24. [5] PUTCHA , N.S. and REDDY , J.N. 1986. A Refind Mixed Shear Flexible F i -nite Element for the Nonlinear Analysis of Laminated Plates. Computer &c Structures, Vol. 22, No. 4, pp. 529-538. [6] REN , J.G. and HINTON, E. 1986. The Finite Element Analysis of Homoge-neous and Laminated Composite Plates Using a Simple Higher Order Theory. Communications in Applied Numerical Methods, Vol. 2, pp. 217-228. [7] T IMOSHENKO, S. and WOINOWSKY-KRIEGER, S. 1959. Theory of Plates and Shells. McGraw-Hill Book Company, New York, U.S.A., pp. 4-15 and pp.367-368. 59 [8] ZIENKIEWICZ, O.C. 1979. The Finite Element Method, 3rd Edition. Tata McGraw-Hill Publishing Company Limited, New Delhi, India, pp. 500-512. 60 Appendix A A Tangential Stiffness Matrix Term The term {a} is required in the derivation of the tangential stiffness matrix [Kt\. Zienkiewicz (1979) presented a procedure to formulate an equivalent expres-sion for this term. Applying this procedure to the analysis, the derivative vector {9} with respect to the natural coordinates £ and f is found to be {H!itS) = {«{P} = lG]{a} where [G] = and {a} is the displacement vector. The nonlinear term in the strain equation expressed in natural coordinates is (A.1) with [A] being Ax t 0 0 & (A.3) Therefore, the term discussed in previous section is now equal to [Bi] = ^[A][G] (AA) Applying a special property described by Zienkiewicz , the term d]£Qy {0} in our beam analysis can now be stated as °~x Ixy 1XV (Ty [G]dV Ax 2 <7x Ixy IXV Oy d{M}T L (o}r dV (A.5) Since ixy and ay are assumed to be O's in the analysis, equation A.5 can be simplified to d[B o T M - AxWv d{M}d{MVdv (A.6) d{a} x2Jv~x d^ d£ This equation can now be substituted into the tangential stiffness matrix equa-tion. 62 Appendix B Material Nonlinearity - Shear Moduli Material nonlinearity occurs because of a nonlinear stress-strain relationship. The discussion in this section will be focussed on the assumption of a nonlinear shear stress-strain relationship for a particular layer. Previous assumptions such as stress continuity and linear shear strain interpolation within each layer remains. B . l General Formulation Consider a general shear stress-strain function r = gfr) (B.l) The elasticity matrix [D] is now dependent upon the strain vector {e}. For example, a commonly used nonlinear stress-strain relationship is given as r = (Po + Pii) (l ~ exp ( ^ 7 ) ) (B.2) B . 2 Piecewise Linear Shear Strain Assumption The shear stress is now a nonlinear function of shear strain for each layer. Recalling the shear strain approximation discussed in chapter 2, the shear strain at 63 the n-th layer is 7 » ( f ) = 7 « . ( ^ ) + 7 » ( ^ ) (B.3) where 7 „ 0 is dependent upon the previous layer's shear strain 7„_i. Shear stress continuity was applied at the interface of any two consecutive layeres. For material nonlinearity and a general stress-strain function, the value 7 n o can now be expressed as ln0 = fn(ln-l) (BA) where f n ( l n - i ) is a nonlinear function of 7„_i. In the previous discussion of linear material properties, f n ( i n - i ) was simply Min-!) = ^ l n - 1 (B.5) However, with a general nonlinear stress-strain relationship, equation B.3 can now be expressed as 7»(f) = / » ( 7 n - l ) ( ^ ) + In ( ^ ) (B.6) The shear distortion u*(£, c) can now be obtained using the same procedure as described in chapter 2. The shear strain in each layer is assumed to be linear and shear stress continuity between layers is imposed. Integrating with respect to the coordinate f , u* at the n-th layer is expressed as -{E^[-r* + / * ( 7 * - i ) ] } fndn-r) + J 64 tNA [LOCAL2 LOCAL l \ l N A [ 4 + 4 + 4j , . . [LOCAL LOCAL2 3\ +fNA[lNA-l) I ~ " + -J (B .7) Expressing this equation in a more compact form, the shear distortion'u*(£, f) can finally be expressed as n NA k=l k=l n NA + J2 T(n, k, f) - E T* (NA> t> LOCAL) (B.8) Jfc=1 *=i where F{n,k,$) h 2 for k = 1,2,... , n — 1 for A; = n F*(NA,k,LOCAL) = 2 for A: = 1,2,... ,n — 1 * A M [LOCAL2 LOCAL ly F(NA, NA, LOCAL) = — — + — + -2 ^ 4 4 4^ for A; = NA 65 and T(n,k,s) = A ( 7 * - i ) | for k = 1,2,...,n— 1 r(»,-,f) = f ^ ' f ^ - ^ l ) for /c = n r'(/\TA,fc, L O C A L ) = /*(7*-i)| for A; = 1 ,2, . . . , n — 1 ^.,»r. x r . r ^ i r , , , .tNA(LOCAL LOCAL2 Z\ T (NA,NA,LOCAL) = f N A d N A - x ) ~ - ( — — + - J for k = NA B .3 Finite Element Formulation The shear distortion u* and shear strain 7 consists of two components: 1. linear component of F and F* and 2. nonlinear component of T and T*. Proceeding as before, the strain equations of du d2w du* l . d t o . j e* = di ~ Zdx* +~dx~ + 2~^~dx~> _ dyS_ lxz — ~~: dz are applied to the displacements. Knowing the shape functions of 66 u = {Nl}T{a} w = {M}T{a} Ik = {N3k}T{a} (B.9) and applying the chain rule, the strain eXn is then expressed as f 2 d{Nl}T ( tn \ 4 d2{M}T E x " ~ \ A x di Vn+ 2 V Ax* di2 ^ 2 „ , , ,d{N3k}T "4 2 , r „ „ ir,d{N3k}T\ f , " 2 a r d7*-i £4 2 ar* d^-i ^ A i a i M de fz\Axdlk^ di ^ • ^ ^ • > <B-l0> The derivative term of 7's respect to i and the partial derivative term of T and r* are then obtained as f i t ! = ^ f ^ W (B.11) a r \ 3/fc(7t-i)«* dlk-i) dln-x 2 for A; = 1,2, . . . ,n — 1 3 / „ ( 7 „ - i ) < » / £ _ £ + 3 57n-i 2 \2 4 4 ; for Jfc = n (B.12) 67 and dlk-i) dik-x 2 for k = l,2t...,NA- 1 [LOCAL LOCAL2 3\ for k = NA (B.13) dfNAJlNA-i) t N A ( 2 , > dlNA-l 2 Equations B . l l to B.13 are now substituted into equation B.10 to give the following bending strain equation eXn = {KX}T{a} + {KXS}T + {a} + {KXG}T{a} (B.14) where tn \ 4 d?{M}T 2 V A x 2 di2 T {KXS} - ^—^—j -{ K X G } r = 1 4 d { M y d { M } r dim*.!? di 2 A x 2 1 1 di di The shear strain has already been determined in the previous section. Further manipulation of equation B.14 gives the shear strain for the n-th layer as 68 = {KXZ}T{a} + {KXZS}T{a} (B.15) where and { K X Z Y = (l±i) {iV3n}> \ 2 / 7„_i Equations B.14 and B.15 are combined to produce the strain vector {e} which is given by {*(»)} = )={B0(n)]{a Zn ) } + [S51(n)]{a} + [S1(n)]{a} (B.16) where [Bo[n)\ = [BS1(n)\ = [Bi(n)} = {KX{n)Y {KXZ(n)}T {KXS{n)Y {KXZS(n)Y ' {KXG(n)Y ' {°y J The matrix [B0] is independent of the displacement vector {a} while the other two matrices [By] and [BSi] are dependent on {a}. The virtual strain equations are now determined in order to produce the virtual work equation. The virtual strain terms from [i?o]{a} and [2?i]{a} have already been determined as 5 {{B0{n)]{a}} = [B0{n)]{Sa} 6 {[BiinMa}} = 2[B1(n)]{6a} = {B2(n)]{6<i} (B.17) Therefore, virtual strain term from [JB<S'1]{a} is now needed to complete the virtual strain equations. Using the chain rule, the term tf{[i?Si]{a}}- is expressed 69 as 6{[BS1{n)}{a}} = 6[BS1[n)]{a} + LBS^n)]^} (B.18) 6{ifX5(n)}r (5{ifXZ,S'(n)}r The term 5{KXS'}r is determined as W + {iOTS(n)}r {iifJ!:Z5(n)}T with £{tfX,S(n)}3 friAx \dlk-i) de d f dT :dlk-i) dak \d~ik-i. / a 2 r \ dt {N3k.1}T{6a} (B.19) which gives (5{KX5(n)}r{a} d f dT ^Ik-i) dak \d% ^ 2 ( d2T de W NA 2 ( a 2 r* ^ A x v a 7 t - 1 i de (B.20) In order to remove the {6a} vector in a later stage, the order of the vectors are rearranged while maintaining the same scalar product. The resulting equation is 6{KXS[n§T{a} f A 2 / a2 r \d{N3^}T T - ^ — { a } { N 3 k - i } 70 £4 2 / d2T* \ d{N3k.1}T , W A r o l T | r c . (B.21) Defining 6{KXS}T{a} = {KXS2?} r {» (B.22) and { K X 5 2 } T = { K X 5 } T + {KXSD}T (B.23) The virtual bending strain can finally be expressed in the finite element form of 6ex = {KX}T{6a} + {KXS2}T{6a} + 2{KXG}T{6a} (B.24) The term 6{KXZS}T is now determined by the same procedure as that used for 6{KXS}T 6{KXZS}T = {NS^f 6 ( / w ^ l ) ) (B.25) where V ln-i J dak ( fnjln-l. V ln-l 1 dfnjln-l) {N3k^}T {6a} (-r^ {N3k.1}T{6a} Therefore ~dfn{ln-l) fnjln-l) 6{KXZS}T{a} = j ^ i - i ^ i N T S ^ } ^ ^ - ) dln-l ln-1 {KXZSD}T{6a} (B.26) 71 and defining {KXZS2}T = {KXZS}T + {KXZSD}1 (B.27) The virtual shear strain is expressed as 6lz*n = {KXZ}T{6a} + {KXZS2}T{6a} Now, the virtual strain vector {6e} can be defined as (B.28) {6e} = | | = [[JBo] + [BS2] + [B2\] {6a} (B.29) where [Bo] = [BS2] = m = {Kxy {KXZ}T {KXS2}T {KXZS2}T 2{KXG}T {0}T Therefore, the column vectors of strain, virtual strain and stress are defined as follow {e} = [Bo + BSl + 5i]{a} {6e} = [B0 + BS2 + B2]{6a} (B.30) w = \Dm and the virtual work equation gives *(o) = [B% + BS* + Bl] [D] [BQ + BSX + Bx] a v ] {a} - {F} = [K\{a}-{F} (B.31) B.4 Newton-Raphson Method If Standard Newton-Raphson method is again used to solve the system of non-linear equations, the tangential stiffness matrix [Kt] must also be determined. From discussion of the Newton-Raphson method in chapter 3, the matrix [Kt] is defined 72 as d{a} d (B.32) d{a} _ d{o) d{e} d{a} ~ d{e} d{a} Substituting this equation into equation B.30 gives Applying chain rule again, the term becomes (B.33) 1 * 1 - L + d [Bl + BS? + Bl) dja} { 0 ) r l d{a} d{c} dV [2?0 + B 5 2 + 5 2 J d { £ } d { a } From previous sections, the term has already been defined as d{s} d{a} = [B0 + BS2 + B2] (B.34) (B.35) The term is the slope of the stress-strain curve at any given strain. This term is then defined as the tangential elasticity matrix [Dt] where [AJ = do, (B.36) The matrix [B0] is independent of the vector {a} thus is equal to zero. The term {a} was already found in chapter 3 as ® W - , [ M I ] (B.37) Finally, d^^yT{<?} has to be determined to complete the [Kt] matrix. From equation B.29, the matrix [ B5 2 ] T is given as 73 [BS2] T _ I T i T (B.38) {KXS2y {KXZS2}T therefore, the terms d{*fay*} and ****g*5a} have to be found. Equation B.23 separated d { * f f a } into two components: d { ™ } S } and • This leads to d {JTXSj} = -^-{KXS} + -^-{KXSD} (B.39) d{a} d{a} d{a} From equation B.20, it was found that d{a} t t l A l W V i / dt " 4 _ 2 _ /_a2r*_ ^ ^ > { / V 3 t _ 1 } r (B.40) d£ However, the term d { K ? { * } D ) T still remains to be found. Differentiating equa-tion B.21 with respect to the vector {a}, the following equation -^{KXSD} = E f / V ( ^ ; - | } r w { A - ) w - , ) r d{a} t = 1 A i 3 7 t . i ^ ^ ( B - 4 1 ) is found. Finally, combining equation B.40 and B.41, the term ^ f f 3 * is expressed as a matrix d {KXS2} = -J^-{KXS} + -J^r{KXSD) = [MXS] (B.42) d{d) d{a} d{a} where 74 [ M X S ] = S s ^ ^ ^ ^ w ^ ^ ^ ' ^ A a ^ l - i di + + " _2_ / d2T \ " 2 / a2r* \ <*e de de de Following the same procedure to find d^K*^y^, the following equations are ob tained. First, the term d^Kd*ayS^ was expressed as {KXZS} = ( ^ — A — A. da dfn{ln-l) _ fnjln-l) ln-1 {N3n.1}{N3n.1}T and then the term d^Kf^yD^ was found to equal 4-{KXZSD) = (±^>l — dfn{ln-\) fn{ln-l) ln-1 { N 3 n _ l } { N 3 n _ l } T + ( l ^ ) [_L. dSfn(ln-l) 2 a/n(7n-l) d7*-i 72_! 57n_! , 2 / n (7„- i ) ' {a}T{N3n^} (B.43) t Adding equations B.43 and B.44 and defining the matrix [MXZS] as [MXZS] = -^-{KXZS}T + ~^-{KXZSD}T d{a} d{a} (B.44) (B.45) 75 where [MXZS] 1 - C \ 1 \dfn{ln-l) fn(ln-l) ( 2 ) [ ln-1 i - <\ r i 0»/ B(7»-i) 2 3/ n(7„-i) + 2/n(7„_i)" 7^-1 ^ r , - l ~3 'n-1 Therefore, d^*jT cr is equal to d[BS2] d{a} -a = oX\MXS) + TXZ\MXZS\ (B.46) and the matrix [Kt] is given by the expression [Kt] = Jv [ax[MXS] + ox[MX] + Txz[MXZS]] [ [5 0 r + BS* + Bl]\Dt][B0 + BS2 + B2]]dV (B.47) It should be noted that the matrices [MX] , [MXS] and [MXZS] are all sym-metric. B .5 Determining the Unknowns From above formulation of [K] and [Kt] matrices, the value of the unknowns fn{ln-i) a n d ^ s derivatives with respect to 7„_i's are required in each iteration. Previous solution of 7's are substituted into the system of nonlinear equations to obtain an improved solution. However, for complicated functions of 9(7), /„(7r»-i) and other unknowns are implicit functions. Thus, the value of the unknowns cannot be directly obtained. The Newton-Raphson method can also be applied to deter-mine these unknowns. Using the stress-strain relationship stated in eqaution B.2 of 76 T = (P 0 + Pa) 1 - exp where PQ, PI and k are known parameters of a layer's material. Applying shear stress continuity between (n-l)-th and n-th layers, the continuity equation .becomes ( i V , + Px^ln-i) ( l - exp [d^^yi = (Po„ + Puln0) (l - exp ( ^ T T n o ) ) (B.48) From this equation, it can be seen that 7 „ 0 (equivalent to / n ( 7 n-i)) cannot be directly obtained. Substituting the material parameters and the previous solution of 7„_i into equation B.48, the left-hand-side of the equation becomes a constant Cn-\. In order to express the continuity equation in a Newton-Raphson format, the constant is move to the right-hand-side of the equation which gives * ( T n o ) = (P0n + Plnlno) (l - exp ( 7 ^ 0 ) ) ~ Cn~l (BA9) Using Newton-Raphson method, the value 7 n o is then determined. Same proce-dure can be applied to determine the other unknowns of a^^"~'^, 9 ^ f f i " - 1 ^ and 77 Appendix C Sample Calculation of Equivalent Elastic Moduli The equivalent elastic moduli (EXtq and EUeil) for the various corrugated layers in a 18 layers (90-26C-42-42-26C-90) beam were calculated. The location of the corrugated layers in the beam was found to have little effect on the computed moduli of the layers. The properties of the 26C corrugation are given as: h = 0.2032 mm. / = 3.8989 mm. / = 3.6068 mm. Ex = 3.1717 x 103MPa Ev = 1.5169 x 10sMPa For a 18 layers beam with 8 layers of corrugation, the distances d's from the midplanes of the corrugated layers to the 2=0 plane are shown in Figure C.l and in Table C.l. Substituting the corrugation's properties into the half wave-length equation (Section 5.2.1) gives t TT2 /3.6068N2] 5 = 3.8989 1 + — L 0„0„ [ 4 V3.8989/ J = 12.1310 mm. 78 2 Z-0-G„ -NA-n-l n + l NL Figure C . l : Distances d's of the 18-layers beam The moment of inertia J in the cross-machine direction is equal to 3.60682(0.2032) 0.6242 mm. 0.81 , 9 c ( 3.6068 V i -+- t.o ^ 3 8 9 8 9 X 2 y For layer #3 the bending rigidties Dx and Dy are respectively found to be: Dx 3.8989 . 0.23023 (3.1717)-12.1310' 12 +3.1717(0.2032) (5.8674) ,3.8989 12.131 7.131 kN-mm 12.131 Dv = 1.5169(0.6242) + 1.5169—^—(0.2032) (5.8674) = 33.963 kN-mm 79 0 Substituting these two values into the equivalent elastic moduli equations, the moduli are found to equal 7.131 JP __ [iSSSi + (3.6O68) (5.8674)2] = 0.05568 kN/mm' 2 - 55.68M Pa _ 33.963 H F - r (3.6068)(5.8674)2] = 0.2652 kN/mm 2 = 265.2MPo Similar calculation for other corrugated layers in the above 18-layers (90-26C-42-42-26C-90) beam gives the following results: Layer # Ex., n (mm.) (MPa) (MPa) 1 13.6906 57.1 265.8 2 9.7790 56.8 265.6 3 5.8674 55.7 265.2 4 1.9558 44.8 260.7 Table C.1: Computed Equivalent Elastic Moduli From table C.1, the distance dn is shown to have little effect on the computed i i equivalent elastic moduli. Therefore, EXeq and Ey<!q are approximated respectively as 55.2 MPa and 262.0 MPa for all corrugated layers. 80 Appendix D Program Listing A listing of the program CUBES is presented in this appendix. The computer language F O R T R A N is used to develop the program. Also, it should be mentioned that a user's manual describing the data file necessary to run the program is avail-able. 81 IMPLICIT REAL*8(A - H,0 - Z) C C THIS PROGRAM PERFORMS A FINITE ELEMENT ANALYSIS C OF A STRUCTURE USING BEAM BENDING ELEMENT WITH SHEAR C INCLUDED. MAXIMUM OF 10 ELEMENTS AND 20 LAMINAE IS ALLOWED. C CUBIC INTERPOLATIONS ARE USED FOR BOTH LONGITUDINAL AND C LATERAL DEFLECTIONS (U AND W). C C DEFINES ALL VARIABLES BEGINNING WITH THE LETTERS A TO H AND C AND 0 TO Z AS DOUBLE PRECISION C C DEFINE VARIABLES C EL CALCULATED ELEMENT LOAD VECTOR FROM PREVIOUS C SOLUTION C SN ELEMENT STIFFNESS MATRIX C ALOAD = STRUCTURE LOAD VECTOR C ICO = ELEMENT CORNER NODE NUMBERS C X GLOBAL COORDINATES OF NODES IN FINITE ELEMENT GRID C IX BOUNDARY CONDITION INTEGER VECTOR, 1 FOR VARIABLE C TO BE RETAINED, 0 FOR ITS ELIMINATION C NO.ALCOPY = MISC. STORAGE C NNODEL= NO. OF NODES PER ELEMENT C NE TOTAL NO. OF ELEMENTS IN STRUCTURE C NNODES= TOTAL NO. OF NODES IN STRUCTURE C NMAT = GROSS NO. OF VARIABLES IN STRUCTURE C NVAR = NO. OF VARIABLES PER NODE C D ELASTIC MODULUS MATRIX C E ELASTIC MODULUS C G SHEAR MODULUS C ESOLTN= ELEMENT SOLUTION VECTOR C PSOLTN= PREVIOUS SOLUTION VECTOR C SOLTN = PRESENT SOLUTION VECTOR C TS0L1 = INCREMENTAL SOLUTION VECTOR C TK = STRUCTURE TANGENT STIFFNESS MATRIX C TN ELEMENT TANGENT STIFFNESS MATRIX C NDIMD = DIMENSION OF THE ELASTIC MODULUS MATRIX, D C NL NUMBER OF LAMINA IN BEAM C NDIMB = OIMENSION OF THE ELEMENTAL MATRICES C TH THICKNESS OF THE LAMINA C Z = Z COORDINATE OF THE MIDPLANE OF THE LAMINA C NPRST= STRESS PRINTOUT CONTROL PARAMETERS C IF NPRST=1, STRESSES ARE PRINTED C IF NPRST=0, NO STRESSES ARE PRINTED C NA = LAMINA WHERE Z=0 AXIS IS LOCATED C LOCAL = LOCAL COORDINATE OF NA-TH LAMINA WHERE C Z=0 AXIS IS LOCATED C DERN1 = VECTOR OF THE 1-ST DERIVATIVE OF U W.R.T. C X-COORDINATE AT EACH INTEGRATION POINT C EDERN3= VECOTR OF THE 1-ST DERIVATIVE OF SHEAR STRAIN C W.R.T. X-COORDINATE AT EACH INTEGRATION POINT C C ALL SYMMETRIC MATRICES ARE STORED AS COLUMN VECTORS C CONTAINING ONLY THE LOWER HALF OF THE MATRICES. C DIMENSION E L ( 4 8 ) , ALOAD(264), ALC0PY(264). AL0ADI(264) DIMENSION IX(264), N0(2), E(20), IC0(1O,2), XX(2) DIMENSION S0LTN(2S4), PSOLTN(264), ES0LTN(48) DIMENSION TN(1176), TK( 12672) 82 DIMENSION T S 0 L K 2 6 4 ) , STRESS(2) COMMON /CONST 1/ D(2,2), G(20), TH(20), DELX, DELY, NDIMD COMMON /C0NST2/ Z(20) COMMON /CONST3/ DERN1(48), EDERN3(48) COMMON /C0NST4/ NL, NDIMB COMMON /C0NST5/ LOCAL, NA COMMON /CONST7/ X(11) DEFINE AS A 1-DIMENSIONAL PROBLEM WITH ONLY ONE BENDING AND ONE SHEAR STRAIN AS UNKNOWNS. NDIMD = 2 SET ITERATION COUNT TO ZERO ITER = O CALL SUBROUTINE TO READ IN ALL FINITE ELEMENT GRID DATA FOR THE PROBLEM CALL LAYOUT(X. ICO, NNODEL, IX, NE, NNODES, NMAT, NVAR) FIND DIMENSION OF ELEMENT ARRAYS FOR THE PROBLEM BEING ANALYSED NDIMB = NNODEL * NVAR INITIALIZE ARRAY TO HAVE ZERO ENTRIES DO 10 J J = 1, NDIMD DO 10 KK = 1, NDIMD ALOADI = INCREMENTAL LOADING VECTOR WHICH STORES THE THE VALUE OF THE LOAD BEING INCREMENTED. D(JJ,KK ) = O.DO DO 20 I = 1, NMAT TS0L1(I) = O.DO ALOAD(I ) = O.DO ALOADI(I) = O.DO SOLTN(I) = O.DO CALL SUBROUTINE TO READ IN THE NUMBER OF LAMINAE. THICKNESS, ELASTIC MODULUS, SHEAR MODULUS, AND Z-COORDINATE OF MID-PLANE OF EACH LAMINA. OTHER PARAMETERS READ IN ARE: THE LOAD, ULTIMATE TENSILE AND COMPRESSIVE CAPACITIES, AND THE STRESS OUTPUT CONTROL VALUE. CALL PROP(E, ALOAD, NVAR. NPRST, NMAT, ALOADI, INCRE) CHECK THE NUMBER OF DEGREES OF FREEDOM PER NODE NDOF = 4 + NL IF (NDOF .NE. NVAR) GO TO 570 PLACE ZEROES AND CONSTANTS INTO VECTOR OF THE 1-ST DERIVATIVE OF SHEAR (EDERN3) W.R.T. X-COORDINATE RESULTED FROM THE NUMBERING SEQUENCE OF THE LAMINAE. 83 DO 30 N2 = 1. NDIMB . 30 EDERN3(N2) = O.DO C DO 70 K = 1, NA IF (K .NE. NA) GO TO 40 EC = 0.5D0 * TH(K) * (L0CAL*O.5D0+0.25D0*L0CAL**2 + 0.25D0) GO TO 60 40 K1 = K + 1 IF (K1 .NE. NA) GO TO 50 EC = 0.5D0 * TH(K) + 0.5D0 * G(K) * TH(K1) * (0.5DO*L0CAL - C 1 25DO*LOCAL**2 + 0.75D0) / G(K1) GO TO 60 50 EC = 0.5D0 * TH(K) + 0.5D0 * G(K) * TH(K1) / G(K1) 60 K4 = K + 4 NEXT = NL + K4 + 4 EDERN3(K4) = -EC * 0.5D0 EDERN3(NEXT) = -EDERN3(K4) 70 CONTINUE C C MAIN PROGRAM C C LHB = LOWER HALF BANDWIDTH OF TANGENTIAL STIFFNESS MATRIX C NA NUMBER OF ELEMENT IN THE VECTOR REPRESENTING THE C TANGENTIAL STIFFNESS MATRIX LHB = NDIMB NA = NMAT * LHB CC CC CHECK FOR INCREMENTAL REQUIREMENTS CC CC IF (INCRE .LE. 1) GO TO 120 C C REMOVE FROM THE LOAD VECTOR THE LOAD BEING INCREMENTED C DO 80 K1 = 1. NMAT 80 ALOAD(K 1 ) = ALOAD(K1) - ALOADI(K1) DO 560 KK = 1, INCRE ITER = 0 C C INCREMENT LOAD C DO 90 MM = 1, NMAT 90 ALOAD(MM) = ALOAD(MM) + ALOADI(MM) / INCRE WRITE (6,100) WRITE (1,100) WRITE (2,100) 100 FORMAT ' ( / , ' WRITE (6,110) KK WRITE (1,110) KK WRITE (2,110) KK 110 FORMAT (/. 1X, 'SOLUTION AT LOAD INCREMENT:', 12) 120 CONTINUE C C INITIALIZE STIFFNESS MATRIX AND A COPY OF THE LOAD VECTOR C TO HAVE ZERO ENTRIES. ALSO ASSIGN EXISTING SOLUTION AS C PREVIOUS SOLUTION AT THE START OF A NEW ITERATION. C DO 130 I = 1, NA 130 T K ( I ) = O.DO 84 DO 140 I = 1. NMAT ALCOPY(I) = 0.00 140 PSOLTN(I) = SOLTN(I) DO 190 I EL = 1, NE DO 180 ILAM = 1. NL DO 160 1 = 1 . NNODEL NO(I) = ICO(IEL.I) DO 150 II « 1. NVAR K = (N0(I ) - 1) * NVAR + 11 L = (I - 1) * NVAR +11 C C IDENTIFY THE PARTICULAR ELEMENT SOLUTION VECTO FROM C THE ENTIRE SOLUTION VECTOR AND IDENTIFY THE X-COORDINATE C OF THE ELEMENT'S NODE. C 150 ESOLTN(L) = PSOLTN(K) 160 XX(I) = X(NO(I)) C C FORM ELASTICITY MATRIX FOR THE ILAM-TH LAYER OF THE C BEAM ELEMENT, D(2 x 2) MATRIX C D( 1 , 1 ) = E(I LAM) D(2,2) = G(ILAM) C C CALCULATE THE ELEMENT'S LENGTH C DELX = XX(2) - XX(1 ) C C CALL SUBROUTINE TO CALCULATE ELEMENT STIFFNESS MATRIX,SN C AND ELEMENT TANGENT STIFFNESS MATRIX.TN USING NUMERICAL C INTEGRATION OF GAUSSIAN QUADRATURE C CALL NSTIFF(I LAM, ESOLTN, EL, TN) C C CALL SUBROUTINE TO INSERT ELEMENT STIFFNESS MATRIX AND C LOAD VECTOR INTO STRUCTURE STIFFNESS MATRIX AND LOAD VECTOR C DO 170 J = 1, NNODEL DO 170 J J = 1, NVAR K = (NO(J) - 1) * NVAR + J J L = ( J - 1) * NVAR + J J C C ADD ALL THE ELEMENTAL CONTRIBUTION TO FORM THE C VECTOR OF LOAD ON THE STRUCTURE CURVE. ALSO PLACE THE C ELEMENT'S TANGENTIAL STIFFNESS MATRIX INTO THE OVERALL C STRUCTURE STIFFNESS MATRIX. C 170 ALCOPY(K) = ALCOPY(K) + EL(L) CALL SETUP(NO, NVAR, TK, TN, NDIMB) 180 CONTINUE 190 CONTINUE C C IMPOSE HOMONGENOUS BOUNDARY CONDITIONS C CALL DISCRD(TK, NMAT, ALCOPY, ALOAD. IX, NDIMB) DO 200 1 = 1 , NMAT C C FIND THE DIFFERENCE BETWEEN THE PRESENT LOAD VALUE AT THE 85 C STRUCTURE CURVE AND THE APPLIED LOAD VALUE. THE CURVE'S C LOAD VALUE IS FOUND BY SUBSTITUTING THE PREVIOUS SOLUTION C VECTOR. THE SYSTEM OF EQUATIONS IS THEN SOLVED USING C CHOLESKY METHOD AND A NEW SOLUTION IS CALCULATED. C 200 T S O L K I ) = ALCOPY(I) - ALOAD (I ) CALL DECOMP(NMAT, NDIMB, TK) CALL SOLV(NMAT, NDIMB, TK. TS0L1) DO 210 I » 1, NMAT 210 SOLTN(I) = PSOLTN(I) - TS0L1(I) C C CALL SUBROUTINE TO CHECK THE ERROR BETWEEN EXISTING C SOLUTION AND PREVIOUS SOLUTION C CALL ERR(PSOLTN, SOLTN, I ERR, NMAT) 220 ITER = 1 + ITER C C WRITE THE DISPLACEMENT VECTOR C IF (ITER .EQ. 1) GO TO 240 IF (I ERR .NE. 0) GO TO 120 WRITE (6,230) 230 FORMAT (/, ' NONLINEAR SOLUTION ) GO TO 260 240 WRITE (6,250) 250 FORMAT (/. ' LINEAR SOLUTION ) 260 WRITE (6,270) ITER 270 FORMAT (/, 'NUMBER OF ITERATION :', 2X. 12) WRITE (6,280) 280 FORMAT (/, 1X, 'SOLUTION VECTOR') C C PRINT OUT THE SOLUTION VECTOR FOR THE LINEAR AND THE C FINAL SOLUTIONS. C NDS = NMAT / NVAR NP = NMAT NDSP = NDS IF (NDS .LE. 6) GO TO 290 NP = 6 * NVAR NDSP = 6 290 WRITE (6,300) (I,1=1,NDSP) 300 FORMAT (/, 2X. 'NODE H' , 2X, 6(12, 12X)) WRITE (6,310) (SOLTN(J),J=1,NP,NVAR) 310 FORMAT (/, 3X, 'W', 3X, 6(2X,G12.4)) WRITE (6,320) (SOLTN(J),J=2,NP,NVAR) 320 FORMAT (/, 3X, 'Wx', 2X, 6(2X.G12.4)) WRITE (6,330) ( SOLTN( <J ) , <J = 3 . NP , NVAR ) 330 FORMAT (/, 3X, 'U'. 3X, 6(2X,G12.4)) WRITE (6,340) (SOLTN(0),J=4,NP,NVAR) 340 FORMAT (/, 3X, 'Ux', 2X, 6(2X,G12.4)) WRITE (6.350) 350 FORMAT (/, 1X, 'SHEAR STRAIN') WRITE (6,360) 360 FORMAT (/. 1X. 'LAMINA H') WRITE (6.370) (SOLTN(J),0=5,NP,NVAR) 370 FORMAT (/, 4X, '1', 2X, 6(2X,G12.4)) IF (NVAR .LT. 6) GO TO 400 DO 380 1 = 6 , NVAR 86 c c c c 1 1 = 1 - 4 380 WRITE (6,390) I I , (S0LTN(J).J=I,NP,NVAR) 390 FORMAT (/, 2X, 13, 2X, 6(2X,G12.4)) PRINT OUT STATEMENTS FOR STRUCTURE WITH MORE THAN SIX NODES. 400 IF (NMAT .LE. NP) GO TO 430 WRITE (6,410) 410 FORMAT (/, ' ') WRITE (6,300) (I.I-7.NDS) NP 1 = NP + 1 WRITE (6,310) (SOLTN(J),d=NP1,NMAT.NVAR) NP2 = NP1 + 1 WRITE (6,320) (SOLTN(0).d=NP2,NMAT,NVAR) NP3 = NP2 + 1 WRITE (6,330) (SOLTN(J).d=NP3.NMAT,NVAR) NP4 = NP3 + 1 WRITE (6,340) (SOLTN(J),J=NP4,NMAT,NVAR) NP5 = NP4 + 1 WRITE (6,350) WRITE (6,360) WRITE (6,370) (SOLTN(J),J=NP5,NMAT,NVAR) IF (NVAR .LT. 6) GO TO 430 NP6 = NP5 + 1 NPVAR = NP + NVAR II = 1 DO 420 I = NP6, NPVAR II = 1 + II 420 WRITE (6,390) I I , (SOLTN(J),d=I,NMAT,NVAR) 430 CONTINUE IF (ITER .EO. 1) GO TO 120 CHECK FOR STRESS PRINTOUT CONTROL 440 450 460 470 480 1 490 500 510 ' ) IF (NPRST .EO. WRITE (6,440) WRITE (1,440) WRITE (2,440) FORMAT (/, 1X, WRITE (6,450) FORMAT (/, 1X, WRITE (6,460) FORMAT (/, 1X, WRITE (6,470) FORMAT (/, 1X, WRITE (1.480) FORMAT (/, 1X, ' WRITE (2,490) FORMAT (/. 1X, DO 550 I EL = 1, WRITE (6,500) WRITE (1,500) WRITE (2.500) FORMAT (/, 1X, WRITE (6,510) FORMAT (/. 1X, 0) GO TO 560 ELEMENT STRESSES') 'ELEMENT STRESSES AT INTEGRATION POINTS ALONG X') 'AT CENTER OF EACH LAMINA') '(BENDING STRESS, SHEAR STRESS)') ' BENDING STRESSES ONLY AT ALL INTEGRATION POINTS SHEAR STRESSES ONLY AT ALL INTEGRATION POINTS') NE I EL I EL I EL 'ELEMENT H', 2X, 13) 'LAMINA H> 2X, 1 ST INTEGRATION POINT', 2X, 87 1 '5 TH INTEGRATION POINT') DO 540 I LAM = 1, NL DO 530 1 = 1 . NNODEL NO(I) = ICO(IEL.I) DO 520 II = 1, NVAR K = (NO(I) - 1) * NVAR + II L = (I - 1) * NVAR + II 520 ESOLTN(L) = SOLTN(K) 530 XX(I) = X(NO(I)) D(1,1) = E(I LAM) D(2.2) = G(ILAM) DELX = XX(2) - XX(1) C C CALL STRESS SUBROUTINE. C CALL FSTRES(ILAM. ESOLTN, I EL, NE. IFAIL) 540 CONTINUE 550 CONTINUE IF ( I F A I L .EO. 1) GO TO 570 560 CONTINUE 570 CONTINUE STOP END C C END OF MAIN PROGRAM C SUBROUTINE LAYOUT(X, ICO, NNODEL, IX, NE, NNODES, NMAT. NVAR) C C THIS SUBROUTINE READS IN ALL FINITE ELEMENT DATA FOR A PROBLEM C C C X ( I ) = X COORDINATES OF I-TH NODE (RETURNED) C ICO(I.d) = U NODE NUMBERS FOR I-TH ELEMENT (RETURNED) C NNODEL = NUMBER OF NODES PER ELEMENT (RETURNED) C IX = BOUNDARY CONDITION CODE VECTOR, =1 IF VARIABLE C IS TO BE RETAINED, =0 IF VARIABLE IS TO BE C ZEROED(RETURNEO) C NE = NUMBER OF ELEMENTS IN TOTAL PROBLEM (RETURNED) C NNODES = TOTAL NUMBER OF NODES IN PROBLEM (RETURNED) C NMAT = GROSS NUMBER OF VARIABLES IN PROBLEM (RETURNED) C NVAR = NO. OF VARIABLES PER NODE (RETURNED) C IMPLICIT REAL*8(A - H,0 - Z) DIMENSION X( 1 1 ) . I C 0 O 0 . 2 ) , IX(264) READ (5,10) NE, NNODES, NVAR, NNODEL 10 FORMAT (415) WRITE (6,20) NE. NNODES, NVAR, NNODEL 20 FORMAT (//. ' NO. OF ELEMENTS '. 15, 5X. 'NO. OF NODES', 15, 5X, 1 'VARIABLES PER NODE', 15, ' NODES PER ELEMENT', 15. /) WRITE (6,30) 30 FORMAT (/, 3X, ' NODE', 7X, ' X-CORD', 6X, ' VARIABLES') DO 60 I = 1, NNODES 12 = NVAR * I 1 1 = 1 2 - NVAR + 1 READ (5,40) X ( I ) . ( I X ( d),J=I 1.12) 40 FORMAT (F12.4, 2412) WRITE (6,50) I, X ( I ) , (IX(d),d=I1,I2) 50 FORMAT (1X, 15. 5X, F12.4. 5X, 2414) 88 60 CONTINUE WRITE (6,70) 70 FORMAT (///, 5X, 'ELEMENT', 5X, 'NODE NUMBERS ') DO 90 I = 1, NE C C READ IN I-TH ELEMENT'S NODE NUMBER C READ (5,80) (IC0(I,d),J=1,NNODEL) 80 FORMAT (1614) 90 WRITE (6,100) I, (ICO(I,J),J=1,NNODEL) 100 FORMAT (5X, 15, 5X, 414, 417) C C FIND THE GROSS NUMBER OF VARIABLES IN STRUCTURE C NMAT = NVAR * NNODES RETURN END C SUBROUTINE PROP(E. ALOAD, NVAR, NPRST, NMAT. ALOADI, INCRE) C C SUBROUTINE TO READ IN THE MATERIAL PROPERTIES, C THICKNESS, AND Z-COORDINATE OF MID-PLANE OF LAMINA C C NL = NUMBER OF LAMINA(RETURNED) C TH = THICKNESS OF LAMINA(RETURNED) C E = ELASTIC MODULUS OF LAMINA(RETURNED) C G = SHEAR MODULUS OF LAMINA(RETURNED) C Z = Z-COORDINATE OF MID-PLANE OF LAMINA(RETURNED) C ALOAD =STRUCTURE LOAD VECTOR (RETURNED) C DELX =ELEMENT LENGTH IN X DIRECTION (RETURNED) C DELY =ELEMENT WIDTH IN Y DIRECTION (RETURNED) C NA =LAMINA NUMBER WHERE Z=0 IS LOCATED (RETURNED) C LOCAL =VALUE IN NATURAL COORDINATE WHERE Z=0 IS C LOCATED(RETURNED) C NPRST =STRESS OUTPUT PARAMETER C IMPLICIT REAL*8(A - H,0 - Z) COMMON /CONST 1/ D(2,2), G(20), TH(20), DELX, DELY, NDIMD COMMON /C0NST2/ Z(20) COMMON /C0NST4/ NL, NDIMB COMMON /C0NST5/ LOCAL, NA COMMON /C0NST6/ ULTCOM, ULTEN COMMON /C0NST7/ X(11) DIMENSION E ( 2 0 ) , AL0AD(264), TL0AD(264), AL0ADI(264) READ (5,10) NL, DELY 10 FORMAT (13, E10.4) WRITE (6,20) NL. DELY 20 FORMAT (/, 1X, 'NUMBER OF LAMINA:', 2X, 13, 6X, 'BEAM WIDTH:', 2X, 1 E10.4) WRITE (6.30) 30 FORMAT (/, 1X, 'LAMINA NUMBER', 2X, 'THICKNESS', 2X, 1 'ELASTIC AND SHEAR MODULI', 2X, 'MID-PLANE Z-COORD') DO 50 I = 1, NL READ (5,40) T H ( I ) , E ( I ) . G ( I ) , Z ( I ) 40 FORMAT (4F12.3) 50 WRITE (6,60) I, T H ( I ) , E ( I ) , G ( I ) , Z ( I ) 60 FORMAT (5X, 13, 8X, F10.4, 3X, F10.1. 2X, F10.1.'7X. F10.4) 89 NPRST •• 1 STRESS OUTPUT IS PRINTED = 0 STRESS OUTPUT IS SUPPRESSED 70 80 90 100 1 10 120 130 140 150 READ (5.70) ULTCOM. FORMAT (2G12.3. 11) WRITE (6,80) FORMAT (/, ' WRITE (6,90) FORMAT (/, ' IF (NPRST .EO. WRITE (6,100) FORMAT (/, '--GO TO 130 WRITE (6,120) FORMAT (/, ' STRESS READ (5.140) NA, LOCAL FORMAT (13, E10.4) WRITE (6.150) NA, LOCAL FORMAT (/, 'Z=0 AXIS IS ULTEN, NPRST ULTCOM ULTIMATE COMPRESSIVE STRENGTH ULTEN ULTIMATE TENSILE STRENGTH = '. 1) GO TO 1 10 NO STRESS OUTPUT = ', G12.4) G12.4) ) IS PRINTED FOR THE FINAL SOLUTION- ') LOCATED IN 1 ' NATURAL COORDINATES ARE:', READ IN CONCENTRATED LOADS ONLY NPLOAD = NO. OF CONCENTRATED LOADS READ (5,160) NPLOAD IF (NPLOAD .EO. 0) GO TO 230 160 FORMAT (13) WRITE (6,170) 170 FORMAT (/. ' WRITE (6.180) 180 FORMAT (/, 1X, 'LOAD # 1 3X, 'LOAD VALUE INCRE = 1 DO 210 11 = 1, NPLOAD READ (5,190) NI, NV, 190 FORMAT (213, E10.4 LAMINA ', 13, 2X, E10.4) CONCENTRATED LOADS 2X, 4X, PLOAD, 12) ' NODE If 'NO. OF INC , 2X, 'DEGREE OF FREEDOM', INCREMENTS') C C C C C C C C C C C C C C C C c c NI = NODE NUMBER WHERE THE CONCENTRATED LOAD IS APPLIED NV = VARIABLE NUMBER OF THE NI-TH NODE WHERE THE CONCENTRATED LOAD IS APPLIED INC = NO. OF LOAD INCREMENTS IF INC = 1 NO INCREMENT IF INC > 1 SAY =3 THEN THE CONCENTRATED LOAD IS DIVIDED INTO 3 EQUAL SEGMENTS AND THE VALUE IS STORED INTO THE PARAMETER INCRP **NOTE** ONLY ONE PARTICULAR LOAD CAN BE INCREMENTED DURING ANY SINGLE RUN OF THE PROGRAM M = NVAR * (NI - 1) + NV STORE DEGREE OF FREEDOM, NUMBER OF INCREMENTS OF THE CONCENTRATED LOAD THAT IS INCREMENTED AND THE VALUE OF THE LOAD INTO SEPARATE VECTOR IF (INC .LE. 1) GO TO 200 INCRE = INC ALOADI(M ) = PLOAD 90 200 ALOAD(M) = PLOAD 210 WRITE (6.220) I I . NI. M, ALOAD(M), INC 220 FORMAT (/, 4X. 12, 4X, 12, 10X. 13, 12X, G12.4, 10X, 12) 230 CONTINUE C C READ IN LATERAL DISTRIBUTED LOAD IN Z DIRECTION WITH C STARTING VALUE OST AT NODEST-TH NODE TO ENDING VALUE OEN C AT NODEN-TH NODE C C NQLOAD = NO. OF DISTRIBUTED LOADS C NODE ST = NODE NUMBER WHERE THE DISTRIBUTED LOAD STARTS C NODEN = NODE NUMBER WHERE THE DISTRIBUTED LOAD ENDS C OST = STARTING LOAD VALUE C OEN = ENDING LOAD VALUE C C DISTRIBUTED LOADS ARE REPLACED BY CONSISTENT LOAD IN C THE PROGRAM C READ (5,240) NOLOAD IF (NOLOAD .EO. O) GO TO 350 240 FORMAT (12) WRITE (6.250) 250 FORMAT (/, 1X. '- — -- DISTRIBUTED LOADS ') DO 340 II = 1, NOLOAD DO 260 L = 1. NMAT 260 TLOAD(L) = O.DO READ (5.270) NODEST. NODEN. OST, OEN, INC 270 FORMAT (213, 2G12.4, 12) WRITE (6,290) WRITE (6,280) I I , NODEST, NODEN, OST, OEN, INC 280 FORMAT (/, 3X. 12, 6X, 2(12,2X). 3X, 2(G12.4,2X), 10X. 12) 290 FORMAT (/, 1X. 'LOAD H', 2X, 'NODE:FROM--TO', 2X, 1 'LOAD VALUE:FR0M--TO', 6X. 'NO. OF INCREMENTS') M2 = (NODEN - 1) * NVAR + 1 M1 = (NODEST - 1) * NVAR + 1 300 LOADEL = NODEN - NODEST ODIFF = (OEN - OST) / LOADEL DO 310 J J = 1, LOADEL 0LOAD1 = OST + ODIFF * ( J J - 1) 0L0AD2 = 0L0AD1 + ODIFF NODE = J J + NODEST - 1 C C TLOAD = TEMPORARY LOAD VECTOR C CALL SUBROUTINE CONSLO TO FIND THE CONSISTENT LOAD C FOR ANY SINGLE ELEMENT WITH LOAD 0L0AD1 AT THE 1-ST C NODE AND LOAD 0L0AD2 AT THE 2-ND NODE. C CALL C0NSL0(0L0A01, 0L0AD2, TLOAD. NODE, NVAR, X) 310 CONTINUE DO 320 1 = 1 . NMAT 320 ALOAD(I) = TLOAD(I) + ALOAD(I) IF (INCRE .GT. 1) GO TO 340 IF (INC .LE. 1) GO TO 340 INCRE = INC C C STORES THE DISTRIBUTED LOAD THAT IS BEING INCREMENTED INTO C A SEPARATE VECTOR 91 c DO 330 J = 1, NMAT 330 ALOADI(J) = TLOAD(J) 340 CONTINUE 350 CONTINUE IF (NPRST .EO. 0) GO TO 380 C C STORES THE FINAL SOLUTION'S BENDING STRESSES AND SHEAR C STRESSES AT ALL INTEGRATION POINTS INTO TWO SEPARATE C FILES. C WRITE (1.360) 360 FORMAT (/. ' STRESS OUTPUT FORMAT H 2 BENDING STRESSES'. ' ONLY') WRITE (2.370) 370 FORMAT (/, ' STRESS OUTPUT FORMAT H 3 SHEAR STRESSES'. ' ONLY') 380 CONTINUE RETURN END C SUBROUTINE CONSLO(QLOAD1, 0L0AD2. TLOAD. NODE. NVAR. X) C C SUBROUTINE TO FIND THE CONSISTENT LOAD TO REPRESENT C THE DISTRIBUTED LOAD C IMPLICIT REAL*8(A - H.O - Z) DIMENSION TL0AD(264), X(11) M1 = (NODE - 1) * NVAR + 1 M2 = M1 + 1 M3 = M1 + NVAR M4 = M3 + 1 NODE 1 = NODE + 1 DELX = X(N0DE1) - X(NODE) TL0AD(M1) = 0.35D0 * DELX * QL0AD1 + 0.15D0 * DELX * 0L0AD2 + 1TLOAD(M1) TL0AD(M2) = (DELX**2) * QL0A01 / 20.D0+(DELX**2) * 0L0AD2 / 30.D0+ 1TLOAD(M2) TL0AD(M3) = 0.15D0 * DELX * 0L0AD1 + 0.35D0 * DELX * QL0AD2 + 1TL0AD(M3) TLOAD(M4) = -(DELX**2) * 0L0AD1 / 30.DO-(DELX**2) * QL0AD2 / 20. 1DO+TL0AD(M4) RETURN END C SUBROUTINE SETUP(NODES. NVAR, A, C, LHB) C C PROGRAM TO SETUP MASTER STIFFNESS MATRIX FROM INDIVIDUAL C FINITE ELEMENT STIFFNESS MATRICES. C C NODES = INTEGER VECTOR CONTAINING ELEMENT NODE NUMBERS C IN ORDER C NVAR = NUMBER OF VARIABLES PER NODE C A = MASTER STIFFNESS MATRIX WHICH MUST BE ZEROED C BEFORE FIRST ENTRY TO SETUP C C = ELEMENT STIFFNESS MATRIX C LHB = LOWER HALF BANDWIDTH INCLUDING DIAGONAL C C IMPLICIT REAL*8(A - H,0 - Z) 92 COMMON /C0NST4/ NL, NDIMB DIMENSION A( 12672) . C(1176). N0DES(2) DO 50 I = 1, NDIMB IM = MOD(I,NVAR) > IF (IM .NE. 0) GO TO 10. NOI = I / NVAR IM = NVAR GO TO 20 10 NOI = (I - IM) / NVAR + 1 20 DO 50 d = 1, I Id = I * (I - 1) / 2 + J dM = MOD(d.NVAR) IF (dM .NE. 0) GO TO 30 NOd = J / NVAR dM = NVAR GO TO 40 30 NOd = (0 - dM) / NVAR + 1 C C FIND POSITIONS (M,N) ROW AND COLUMN OF MATRIX A. THE C STRUCTURE'S TANGENTIAL STIFFNESS MATRIX, AND PLACES C INTO THESE POSITIONS THE CORRESPONDING VALUE FROM C, THE C ELEMENT'S TANGENTIAL STIFFNESS MATRIX. 40 M = (NODES(NOI) - 1) * NVAR + IM N = (NODES(NOd) - 1) * NVAR + dM MN = (LHB - 1) * (N - 1) + M 50 A(MN) = A(MN) + C ( I d ) RETURN END C SUBROUTINE DISCRD(TK, NMAT. ALCOPY, ALOAD, IX. LHB) C C THIS SUBROUTINE APPLIES HOMOGENEOUS BOUNDARY CONDITIONS BY C PLACING ZEROS ON OFF DIAGONAL TERMS AND ONE ON THE DIAGONAL C TERM OF THE STIFFNESS MATRICES. ALSO THE LOAD VECTOR TERM C IS REPLACED BY THE HOMOGENEOUR BOUNDARY VALUE OF ZERO. C C TK =TANGENT STIFFNESS MATRIX (RETURNED) C NMAT =GROSS SIZE OF PROBLEM (GIVEN) C ALOAD =GIVEN LOAD VECTOR (RETURNED) C ALCOPY CALCULATED LOAD VECTOR (RETURNED) C IX =BOUNDARY CONDITION CODE VECTOR (GIVEN) C LHB =HALF BANDWIDTH INCLUDING THE DIAGONAL TERM (GIVEN) C IMPLICIT REAL*8(A - H,0 - Z) DIMENSION IX(264), AL0AD(264), ALC0PY(264). TK(12672) MULTI = LHB - 1 DO 50 L = 1, NMAT IF ( I X ( L ) .NE. O) GO TO 50 Id = MULTI * L + L - MULTI TK(Id) = 1.D0 ALCOPY(L) = O.DO ALOAD(L) = O.DO DO 10 N = 1, MULTI IdC = N + Id 10 TK(IdC) = O.DO IF (L .LT. LHB) GO TO 30 DO 20 N1 = 1, MULTI IdR = Id - MULTI * N1 93 20 TK(IJR) = O.DO GO TO 50 30 LL = L - 1 DO 40 N2 = 1, LL IJR = IJ - MULTI * N2 40 TK(IJR) = O.DO 50 CONTINUE RETURN END C SUBROUTINE ERR(PSOLTN, SOLTN, I ERR, NMAT) C C SUBROUTINE TO CHECK THE NONLINEAR SOLUTION HAS C CONVERGED TO DESIRED ACCURACY IN PERCENT C C PSOLTN =PREVIOUS SOLUTION (GIVEN) C SOLTN =PRESENT SOLUTION (GIVEN) C I ERR =ERROR INDICATOR : O ERROR COMPARISON OK (RETURNED) C 1 ERROR COMPARISON FAILED C NMAT =GROSS NUMBER OF EQUATIONS WITH BOUNDARY CONDITIONS C ELIMINATED C ALPER ^ALLOWABLE ERROR PERCENTAGE C ABPER =ABSOLUTE VALUE OF CALCULATED ERROR PERCENTAGE C IMPLICIT REAL*8(A - H.O - Z) DIMENSION PS0LTN(264), S0LTN(264) I = 0 I ERR = O C C SET THE ALLOWABLE ERROR PERCENTAGE, ALPER C ALPER = 0.01D0 10 I = I + 1 C C CHECK THAT SOLUTION IS NOT EQUAL TO ZERO OR LESS THAN C 10E-6 IN ABSOLUTE VALUE. IF THE ABOVE IS TRUE THEN THE C ERROR CALCULATION IS IGNORED. C IF (DABS(SOLTN(I ) ) . L E . 0.000001D0) GO TO 20 DIFF = PSOLTN(I) - SOLTN(I) PER = DABS(DIFF) / DABS(SOLTN(I)) IF (PER .GT. ALPER) GO TO 30 20 CONTINUE IF (I .EQ. NMAT) GO TO 40 GO TO 10 30 I ERR = 1 40 CONTINUE RETURN END C SUBROUTINE STRKX(II, vJJ, I LAM, VALUE, AKX, VALUE3) C C SUBROUTINE TO GENERATE THE VECTOR (AKX) r* Z II 'LOCATION OF INTEGRATION POINTS ALONG Z DIRECTION : UJ 'LOCATION OF INTEGRATION POINTS ALONG X DIRECTION (GIVEN) 3 I LAM =LAMINA NUMBER (GIVEN) ; VALUE=GAUSSIAN INTEGRATION POINTS (GIVEN) 94 C AKX =LINEAR STRAIN VECTOR IN X DIRECTION (RETURNED) C IMPLICIT REAL*8(A - H.O - Z) COMMON /CONST 1/ D(2,2), G(20), TH(20), DELX, DELY. NDIMD COMMON /CONST2/ Z(20) COMMON /C0NST3/ DERN1(48), EDERN3(48) COMMON /C0NST4/ NL, NDIMB DIMENSION AKX(48) , VALUE(5), DDERM(48), FDERN3(48), TEMP(1) DIMENSION VALUE3(3) CALL BSTRAM(JJ, VALUE, DELX. DDERM) CALL BSTRAF(II, VALUE3, I LAM, FDERN3) CALL BSTRAU(Jd, VALUE, DELX, DERN1) C C LINEAR PORTION OF BENDING STRAIN OF THE C STIFFNESS MATRICES C C DERN1 = VECTOR OF THE VALUE dU/dX C DDERM = VECTOR OF THE VALUE d2W/dX2 C FDERN3= VECTOR OF THE X-DEPENDENT VALUE OF dU*/dX C EDERN3= VECTOR OF THE CONSTANT VALUE OF dU*/dX C DO 10 K = 1, NDIMB TEMPO) = (Z(ILAM) + 0.5D0*TH(I LAM)*VALUE3(11)) * DDERM(K) * 2. 1 DO / DELX 10 AKX(K) = 2.DO * (DERN1 (K) - TEMPO) + FDERN3(K) - EDERN3(K)) / 1DELX RETURN END C SUBROUTINE BSTRAM(JJ, VALUE, DELX, DDERM) C C SUBROUTINE TO GENERATE THE 2ND DERIVATIVE C OF THE VECTOR (M) C C dd.VALUE,DELX =SEE PREVIOUS EXPLANATIONS (GIVEN) C ODERM =VECTOR CONTAINING THE 2 ND DERIVATIVE C OF THE FUNCTION RELATING THE TRANSVERSE C DISPLACEMENT, W, AND THE NODAL DISPLACEMENTS C (M) (RETURNED) C IMPLICIT REAL*8(A - H.O - Z) COMMON /C0NST4/ NL, NDIMB DIMENSION VALUE(5) , DDERM(48) DO 10 L = 1, NDIMB 10 DDERM(L) = O.DO DDERMO) = VALUE (dd) * 1 . 5D0 DDERM(2) = (0.75D0*VALUE(dd) - 0.25D0) * DELX L1 = 5 + NL L2 = 6 + NL DDERM( L 1 ) = -DDERMO) DDERM(L2) = DDERM(2) + 0.5D0 * DELX RETURN END C SUBROUTINE BSTRAF(II, VALUES, I LAM, FDERN3) C C SUBROUTINE TO GENERATE THE DERIVATIVE OF THE C VECTOR {N3}XF 95 C 11,VALUE3,I LAM 'SEE PREVIOUS EXPLANATIONS (GIVEN) C FDERN3 'VECTOR CONTAINING THE 1 ST DERIVATIVE OF THE C FUNCTION RELATING THE X-DISPLACEMENT DUE TO C SHEAR DISTORTION OF CROSS-SECTION ANO THE C NODAL DISPLACEMENT (N3) (GIVEN) C IMPLICIT REAL*8(A - H.O - Z) COMMON /CONST 1/ D(2.2), G(20). TH(20). DELX. DELY. NDIMD COMMON /C0NST4/ NL, NDIMB DIMENSION FDERN3(48), VALUE3(3) DO 10 KK = 1, NDIMB 10 FDERN3(KK) = O.DO DO 50 K = 1 . I LAM IF (K .NE. I LAM) GO TO 20 F = 0.5D0 * TH(K) * (VALUES(II)*0.5D0+0.25D0*VALUE3(II)**2 + O. 1 25D0) GO TO 40 20 KI = K + 1 IF (K1 .NE. I LAM) GO TO 30 F = 0.5D0 * TH(K) + 0.5D0 * G(K) * TH(K1) * (0.5D0*VALUE3( II) -1 0.25D0*VALUE3(II)**2 + 0.75DO) / G(K1) GO TO 40 30 F = O.SDO * TH(K) + 0.5D0 * G(K) * TH(K1) / G(K1) 40 K4 = K + 4 NEXT = NL + K4 + 4 FDERN3(K4) = -F * 0.5D0 FDERN3(NEXT) = -FDERN3(K4) 50 CONTINUE 60 CONTINUE RETURN END C SUBROUTINE STRKXZ(II, dd, I LAM. VALUE3, AKXZ. VALUE) C C SUBROUTINE TO GENERATE VECTOR (AKXZ) C C LINEAR PORTION OF SHEAR STRAIN OF THE C STIFFNESS MATRICES C C II.dd.ILAM,VALUE3 =SEE PREVIOUS EXPLANATIONS (GIVEN) C AKXZ 'LINEAR SHEAR STRAIN VECTOR IN C X-Z PLANE(RETURNED) C IMPLICIT REAL*8(A - H.O - Z) COMMON /CONST 1/ D(2,2). G(20), TH(20). DELX, DELY, NDIMD COMMON /C0NST4/ NL, NDIMB DIMENSION AKXZ(48), S1N3(48), S2N3(48), VALUE3(3), VALUE(5) CALL SSTRA(I LAM, VALUE, G. dd, S1N3) IF (I LAM .GT. 1) GO TO 20 DO 10 M = 1, NDIMB 10 S2N3(M) = O.DO GO TO 30 20 CONTINUE LOWER = I LAM - 1 CALL SSTRA(LOWER, VALUE, G, dd, S2N3) 30 CONTINUE DO 40 MM = 1, NDIMB 96 40 AKXZ(MM) = ((1 + VALUE3(II))*S1N3(MM) + (1 - VALUE3(II))*S2N3(MM)) 1 / (2*G(ILAM)) RETURN END C SUBROUTINE SSTRA(NLAM. VALUE. G, J J . SN3) C C SUBROUTINE OF STRKXZ C GENERATE THE VECTOR <N3}XG(NLAM) C . C NLAM =LAMINA NUMBER (GIVEN) C dd,VALUE=SEE PREVIOUS EXPLANATIONS (GIVEN) C G =SHEAR MODULUS (GIVEN) C SN3 =VECTOR OF {N3> TIMES SHEAR MODULUS (RETURNED) C IMPLICIT REAL*8(A - H,0 - Z) COMMON /C0NST4/ NL. NDIMB DIMENSION SN3(48), VALUE(5). G(20) NL1 = NLAM + 4 NL2 = NL1 + NL + 4 DO 10 K = 1, NDIMB 10 SN3(K) = O.DO SN3(NL1) = (1.DO-VALUE(dd)) * G(NLAM) * 0.5D0 SN3(NL2) = (1.D0+VALUE(dd)) * G(NLAM) * 0.5D0 RETURN END SUBROUTINE NSTIFF(I LAM, ESOLTN, EL, TN) C C SUBROUTINE TO CALCULATE THE LINEAR COMPONENT OF C THE STIFFNESS MATRIX. [BO]T[D][BO], USING A C 5 POINTS GAUSSIAN INTEGRATION FOR AN ELEMENT C C I LAM =LAMINA NUMBER (GIVEN) C ESOLTN=PREVIOUS SOLUTION (GIVEN) C TN =NONLINEAR ELEMENT STIFFNESS MATRIX (RETURNED) C AKX,AKXZ.B0.B1,B2=TEMP0RARY STORAGE(GENERATED) C TEMP2,TEMP3,TEMP4=TEMP0RARY STORAGE(GENERATED) C IMPLICIT REAL*8(A - H,0 - Z) COMMON /CONST 1/ D(2,2), G(20). TH(20), DELX, DELY, NDIMD COMMON /C0NST2/ Z(20) COMMON /C0NST4/ NL, NDIMB DIMENSION TN(117G), TE1(2,48), TE2(2,48), WEIGHT(5), ES0LTN(48) DIMENSION VALUE(5), AKX(48), AKXZ(48), AMX(48,48), B2R(48) DIMENSION TEMP 1 (48,48), TEMP2(48,48), EL(48), SN(48,48) DIMENSION B0(2,48), B1(2,48), B2(2,48), BOK2.48). B02(2.48) DIMENSION STRAIN(2), STRESS(2), WEIGH3(3), VALUE3(3) DATA VALUE /-0.90G1798459, -0.5384693101, 0.000000000000000, 1 0.5384693 101. 0.906 1798459/ DATA WEIGHT /O.2369268850, 0.4786286705, 0.5688888889, 1 0.4786286705, 0.2369268850/ DATA VALUE3 /-0.7745966692, 0.000000000, 0.7745966692/ DATA WEIGH3 /O.5555555555. 0.888888889. 0.5555555555/ NT = NDIMB * (NDIMB + 1) / 2 C C TN ELEMENT TANGENTIAL STIFFNESS MATRIX C (SYMMETRIC) C SN ELEMENT STRUCTURAL STIFFNESS MATRIX 97 C (NONSYMMETRIC) C DO 10 13 « 1. NT 10 TN(13) = O.DO DO 20 I 1 = 1. NDIMB DO 20 12 = 1, NDIMB 20 SN(I 1,12) = O.DO DO 80 I = 1, 3 00 80 d = 1 , 5 C C I IS FOR PHI(IN Z) VALUES, AND d IS FOR PSI(IN X) VALUES C CALL STRKX(I, d, I LAM, VALUE, AKX. VALUE3) CALL STRKXZ(I, d, I LAM, VALUES, AKXZ. VALUE) CALL NBSTMX(d, VALUE, DELX, AMX) C C FORM THE [BJ'S MATRICES FROM THE VECTORS AKX. AKXZ, AMX, C AND THE SOLUTION VECTOR C DO 40 dd = 1. NDIMB TEM1 = O.DO DO 30 II = 1. NDIMB TEM1 = ESOLTN(II) * AMX(II.dd) + TEM1 30 CONTINUE B2R(dd) = TEM1 40 CONTINUE 00 50 L = 1, NDIMB B1( 1 ,L) = 0.5D0 * B2R(L) B 1 (2,L ) = O.DO B2( 1 ,L) = B2R(L) B2(2,L) = O.DO B0(1,L) = AKX(L) B0(2,L) = AKXZ(L) B01(1,L) = B0(1.L) + B1(1,L) B01(2,L) = B0(2.L) + B1(2,L) B02(1,L) = B0(1,L) + B2(1,L) 50 B02(2.L) = B0(2,L) + B2(2,L) C C FORM THE ELEMENT'S TANGENTIAL STIFFNESS MATRIX, TN C CALL DGMULT(D, B02, TE 1 , NDIMD, NDIMD, NDIMB, NDIMD, NDIMD, 1 NDIMD) CALL DGPR0D(B02, TE1, TEMPI, NDIMB. NDIMD, NDIMB. 2. 2, 48, 1 1 O. 1) CALL DGMATV(B01, ESOLTN, STRAIN, NDIMD, NDIMB, NDIMD) CALL DGMATV(D. STRAIN. STRESS. NDIMD, NDIMD, NDIMD) DO GO M = 1, NDIMB DO SO N = 1. M M N = M * ( M - 1 ) / 2 + N 60 TN(MN) = WEIGH3(I) * WEIGHT(d) * DELX * 0.25DO * TH(ILAM) * 1 DELY * (TEMPI(M,N) + STRESS(1 ) *AMX(M,N)) + TN(MN) C C FORM THE ELEMENT'S STRUCTURE STIFFNESS MATRIX, SN C CALL DGMULT(D, B01 , TE2, NDIMD, NDIMD, NDIMB, NDIMD, NDIMD, 1 NDIMD) CALL DGPR0D(B02, TE2, TEMP2, NDIMB. NDIMD, NDIMB, 2, 2, 48, 1 1 0. 2) 98 DO 70 M = 1, NDIMB DO 70 N = 1, NDIMB 70 SN(M,N) = WEIGH3(I) * WEIGHT(d) * DELX * 0.25D0 * TH(ILAM) * 1 DELY * TEMP2(M,N) + SN(M.N) 80 CONTINUE C C FORM THE LOAD CURVE VALUE VECTOR, EL WHERE EL IS THE C VECTOR OF LOAD VALUES ON THE CURVE WITH THE PREVIOUS C SOLUTION SUBSTITUTED INTO THE SYSTEM OF EQUATIONS. C DO 100 II = 1, NDIMB TEM = O.DO DO 90 dd = 1. NDIMB 90 TEM = SN(II.dd) * ESOLTN(dd) + TEM E L ( I I ) = TEM 100 CONTINUE RETURN END C SUBROUTINE NBST'MX (d J , VALUE, DELX, AMX ) C C SUBROUTINE TO GENERATE THE NON-LINEAR PART C OF THE STIFFNESS MATRIX. [MX] C C dd,VALUE.DELX =SEE PREVIOUS EXPLANATIONS (GIVEN) C AMX =NONLINEAR STRAIN MATRIX OF THE 1 ST DERIVATIVE C OF {M> TIMES ITS TRANSPOSE (RETURNED) C IMPLICIT REAL*8(A - H.O - Z) COMMON /C0NST4/ NL, NDIMB DIMENSION VALUE(5), AMX(48,48), DM(48), VALUE3(3) DO 10 K = 1, NDIMB 10 DM(K) = O.DO DM(1) = -1.5D0 - (VALUE(dd) + 1.D0) + 0.75DO • (VALUE(dd) + 1.00) 1** 2 DM(2) = (0.5D0-(VALUE(dd) + 1.D0) + 0.375DO*(VALUE(dd) + 1.D0)**2) 1 * DELX K1 = NL + 5 K2 = NL + 6 DM(K1) = -DM(1) DM(K2) = (-0.5D0*(VALUE(dd) + 1.D0) + 0.375D0*(VALUE(dd) + 1.D0)** 12) * DELX DO 20 I = 1, NDIMB DO 20 d = 1 , NDIMB 20 AMX(I.d) = DM(I) * DM(d) * 4.DO / (DELX**2) RETURN END C SUBROUTINE DECOMP(N, LHB, A) C C SUBROUTINE TO DECOMPOSE A MATRIX IN A SYSTEM OF EQUATIONS C USING CHOLESKY METHOD FOR BANDED, SYMMETRIC. POSITIVE C DEFINITE MATRIX TO SOLVE THE SYSTEM. C IMPLICIT REAL*8(A - H.O - Z) DIMENSION A(12672) C C A IS STORED COLUMNWISE 99 c KB = LHB - 1 C C DECOMPOSITION C A( 1 ) = DSORT(A(1) ) IF (N .EO. 1) RETURN DO 10 I = 2 . LHB 10 A ( I ) = A ( I ) / A(1) DO 60 J = 2, N J1 = d - 1 IdD = LHB * J - KB SUM = A(IdD) KO = 1 IF (d .GT. LHB) KO = d - KB DO 20 K = KO. d l JK = KB * K + d - KB 20 SUM = SUM - A(dK) * A(dK) A(IJD) = DSORT(SUM) DO 50 I = 1. KB II = d + I KO = 1 IF ( I I .GT. LHB) KO = II - KB SUM = A(IdD + I ) IF (I .EO. KB) GO TO 40 DO 30 K = KO. d1 dK = KB * K + d - KE IK = KB * K + 11 - KB 30 SUM = SUM - A(IK) * A(dK) 40 A( IdD + I) = SUM / A(IdD) 50 CONTINUE 60 CONTINUE RETURN END C SUBROUTINE SOLV(N, LHB, A, B) C C SUBROUTINE TO SOLVE THE SYSTEM OF EQUATIONS USING CHOLESKY C METHOD AFTER THE MATRIX HAS BEEN DECOMPOSED BY DECOMP C SUBROUTINE CAN SOLVE FOR DIFFERENT VALUES OF B WITHOUT C DECOMPOSING THE MATRIX A REPEATEDLY. C IMPLICIT REAL*8(A - H.O - Z) DIMENSION A(12672), B(264) C C FORWARD SUBSTITUTION C KB = LHB - 1 B(1) = B(1) / A(1) IF (N .EO. 1) GO TO 30 DO 20 I = 2, N 1 1 = 1 - 1 KO = 1 IF (I .GT. LHB) KO = I - KB SUM = B ( I ) 1 1 = LHB * I - KB DO 10 K = KO, 11 IK = KB * K + I - KB 100 10 SUM = SUM - A(IK) * B(K) B ( I ) = SUM / A ( I I ) 20 CONTINUE C C BACKWARD SUBSTITUTION C 30 N1 = N - 1 LB = LHB * N - KB B(N) = B(N) / A(LB) IF (N .EO. 1) RETURN DO 50 I = 1, Nl I 1 = N - I + 1 NI = N - I KO = N IF (I .GT. KB) KO = NI + KB SUM = B(NI) II = LHB * NI - KB DO 40 K = 11, KO IK = KB * NI + K - KB 40 SUM = SUM - A ( I K ) * B ( K ) B(NI) = SUM / A ( I I ) 50 CONTINUE RETURN END C SUBROUTINE BSTRAU(JJ. VALUE, DELX, DERN1) C C SUBROUTINE TO GENERATE THE 1ST DERIVATIVE OF THE VECTOR {U} C C JJ.VALUE,DELX =SEE PREVIOUS EXPLANATION C DERN1 =VECTOR CONTAINING THE 1ST DERIVATIVE OF C THE AXIAL DISPLACEMENT, U, AND THE NODAL C DISPLACEMENT VECTOR. (RETURNED) C IMPLICIT REAL*8(A - H.O - Z) COMMON /C0NST4/ NL, NDIMB DIMENSION VALUE(5), DERN1(48) DO 10 L = 1, NDIMB 10 DERNKL) = O.DO DERN1(3) = -1.5D0 * (VALUE(JJ) + 1.D0) + 0.75D0 * (VALUE(JJ) + 1. 1D0) ** 2 DERN1(4) = (O.5D0-(VALUE(JJ) + 1.D0) + 0.375D0*(VALUE(JJ) + 1.00)* 1*2) * DELX L1 = NL + 7 L2 = NL + 8 DERN1(L 1 ) = -DERN1(3) DERN1(L2) = (-0.5D0*(VALUE(Jd) + 1.D0) + O.375D0*(VALUE(JJ) + 1. 1D0)**2) * DELX RETURN END C SUBROUTINE FSTRES(ILAM, ESOLTN, IEL, NE, IFAIL) C C SUBROUTINE TO CALCULATE THE STRESS OF AN ELEMENT C USING A 5 POINTS GAUSSIAN INTEGRATION IN X C AND A 3 POINTS INTEGRATION IN Z C C ILAM =LAMINA NUMBER (GIVEN) 101 c c c c c c c c ESOLTN=PREVIOUS SOLUTION (GIVEN) AKX,AKXZ,BO,B1,B2=TEMPORARY STORAGE(GENERATED) TEMP2,TEMP3,TEMP4=TEMP0RARY STORAGE(GENERATED) IMPLICIT REAL*8(A - H.O - Z) COMMON /CONST 1/ D(2,2), G(20), TH(20), DELX, DELY, NDIMD COMMON /CONST2/ Z(20) COMMON /CONST4/ NL, NDIMB COMMON /C0NST6/ ULTCOM, ULTEN DIMENSION ESOLTN(48) DIMENSION VALUE(5), AKX(48), AKXZ(48), AMX(48,48), B2R(48) DIMENSION BO(2,48), B 1 (2,48 ) , B01(2,48) DIMENSION STRAIN(2), STRESS(2). BENDS(5), SHEARS(5), VALUE3(3) DATA VALUE /-0.9061798459, -0.5384693101. 0.OOOOOOOOOOOOOOO. 1 0.5384693101, O.9061798459/ DATA VALUE3 /-0.7745966692, 0.0000000000. 0.7745966692/ IF (I EL .NE. 1) GO TO 10 TENMAX = O.DO COMMAX = O.DO 10 WRITE (2,20) I LAM WRITE (1,20) I LAM 20 FORMAT (/, ' LAMINA H '. 12) DO 140 I = 1, 3 DO 90 J = 1, 5 I IS FOR PHI(IN Z) VALUES, AND J IS FOR PSI(IN X) VALUES CALL STRKX(I, J , I LAM, VALUE, AKX, VALUE3) CALL STRKXZ(I. J . I LAM, VALUES, AKXZ. VALUE) CALL NBSTMX(J, VALUE, DELX, AMX) DO 40 UJ = 1, NDIMB TEM1 = O.DO DO 30 II = 1, NDIMB TEM1 = ESOLTN(II) CONTINUE B2R(JJ) = TEM1 * AMX(II.JJ) + TEM1 30 40 B2R(L) NDIMB 0.5D0 * O.DO AKX(L) AKXZ(L) BO(1.L) + B1(1,L) CONTINUE DO 50 L = 1 B1(1,L) = B1 (2 , L) = B0(1,L) = B0(2,L) = B01(1,L) = 50 B01(2,L) = B0(2,L) + B1(2,L) CALL DGMATV(BO1, ESOLTN, STRAIN, NDIMD, NDIMB, NDIMD) CALL DGMATV(D, STRAIN, STRESS, NDIMD, NDIMD, NDIMD) CALCULATE THE BENDING AND SHEAR STRESSES AT THE PARTICULAR INTEGRATION POINT BENDS(J) = STRESS(1) SHEARS(J) = STRESS(2) FIND THE MAXIMUM COMPRESSIVE AND TENSILE STRESSES FOR THE BEAM AND RECORD THEIR LOCATIONS IF (BENDS(J) .LE. O.DO) GO TO 70 IF (BENDS(J) .LE. TENMAX) GO TO 60 102 TENMAX = BENDS(d) MAXTEL = 1 EL MAXTLA = ILAM MAXTI = I MAXTd = d 60 IF (BENDS(d) .LT. ULTEN) GO TO 90 C C SET ANY TENSILE STRESS THAT EXCEEDS THE ALLOWABLE TENSILE C STRESS AS EQUAL TO 1.00E+30 C BENDS(d) = 1.00E+30 IFAIL = 1 GO TO 90 70 IF (BENDS(d) .GE. COMMAX) GO TO 80 COMMAX = BENDS(d) MAXCEL = I EL MAXCLA = ILAM MAXCI = I MAXCd = d 80 IF (DABS(BENDS(d) ) .LT. ULTCOM) GO TO 90 C C SET ANY COMPRESSIVE STRESS THAT EXCEEDS THE ALLOWABLE C COMPRESSIVE STRESS AS EQUAL TO -1.00E+30 C BENDS(d) = -1..00E + 30 IFAIL = 1 90 CONTINUE C C PRINT OUT THE STRESS IF REQUESTED C IF (NL .EQ. 1) GO TO 100 IF (I .NE. 2) GO TO 120 100 WRITE (6.110) ILAM, BENDS(1), SHEARS(1), BENDS(5), SHEARS(5) 110 FORMAT (/, 13, 6X, 2('(',E9.3, ', ' , E9.3, ' ) ' ,3X)) 120 WRITE (1,130) (BENDS(II),II=1.5) 130 FORMAT (/, 2X, 5(G12.4,1X)) WRITE (2.130) (SHEARS(II),11=1,5) 140 CONTINUE IF ( ( I EL .NE. NE) .OR. (ILAM .NE. NL)) GO TO 220 WRITE (6,150) MAXCEL 150 FORMAT (/. ' MAXIMUM COMPRESSION OCCURS AT ELEMENT H '. 12) WRITE (6,160) MAXCLA 160 FORMAT (1X, 'LAMINA H ', 12) WRITE (G.170) MAXCI, MAXCd 170 FORMAT (1X, 'INTEGRATION POINTS I AND d OF ', 12, 2X, 'AND ', WRITE (6,180) COMMAX 180 FORMAT (1X, 'COMPRESSION = ', G12.4) WRITE (6,190) MAXTEL 190 FORMAT (/. ' MAXIMUM TENSION OCCURS AT ELEMENT H '. 12) WRITE (6,160) MAXTLA WRITE (6,170) MAXTI, MAXTd WRITE (6,200) TENMAX 200 FORMAT (1X, 'TENSION = '. G12.4) IF (IFAIL .NE. 1) GO TO 220 WRITE (6.210) 210 FORMAT (/. '********* THE BEAM HAS FAILED. *******••') 220 CONTINUE RETURN END 103
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multilayer beam analysis including shear and geometric...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multilayer beam analysis including shear and geometric nonlinear effects Kwan, Herman Ho Ming 1987
pdf
Page Metadata
Item Metadata
Title | Multilayer beam analysis including shear and geometric nonlinear effects |
Creator |
Kwan, Herman Ho Ming |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | This thesis presents an analysis and experimental verification for a multilayer beam in bending. The formulation of the theoretical analysis includes the combined effect of shear and geometric nonlinearity. From this formulation, a finite element program (CUBES) is developed. The experimental tests were done on multilayer, corrugated paper beams. Failure deflections and loads are thus obtained. The experimental results are reasonably predicted by the numerical results. Based upon this comparison, a maximum compressive stress is determined for the tested beam. Finally, design curves for the tested beam are drawn using the determined maximum compressive stress and the finite element program. |
Subject |
Girders -- Testing Shear (Mechanics) |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062836 |
URI | http://hdl.handle.net/2429/26711 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A7 K87_5.pdf [ 7.1MB ]
- Metadata
- JSON: 831-1.0062836.json
- JSON-LD: 831-1.0062836-ld.json
- RDF/XML (Pretty): 831-1.0062836-rdf.xml
- RDF/JSON: 831-1.0062836-rdf.json
- Turtle: 831-1.0062836-turtle.txt
- N-Triples: 831-1.0062836-rdf-ntriples.txt
- Original Record: 831-1.0062836-source.json
- Full Text
- 831-1.0062836-fulltext.txt
- Citation
- 831-1.0062836.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062836/manifest