. .1 33.21.91.75... . . 72313"I1i. iii. .5— : V . , _‘ L5, 3311' f 3 It. "I" I. “L t 1 nu . DI ‘ {k1 11'3_ — ’ —:— filly". .L ' .' L. I #123ij I“ 1* s .664 1.2 IE £13 2' 3 'I 3” r- A.’ , A 1 {:5 :11;‘3 M 5139' a! 11‘. 3'3... 1L 1} 1 .1 II J g. I I: . 2? 11’1" . _§ 33313 £12. ‘( I7 I: , l 1‘, 4 | ‘ G- I l at“, .. ' Kij : ' AG ‘1'- v 1 W323! 31:31:; .3; o. I j)? 33*?111431 ‘JHILU MM ‘ p.411? :3}. oh " 2 ‘ ”3W! 4‘|V1‘ 1 I ’8 ‘|L"4Y:?‘\n‘¢t! .J_ “,3“. 3'33 n‘:- ..~..’1%r‘|7,rrc.€k '1....‘.|\‘(45§{JJ) 1121.213}?! H.355 K' . .1 .1 1‘1 1 "’1’ 33 i'uL -'!‘1."3"1 ‘1’ 1."? 1, 1 .1013)" f‘.-¢ fill-{$.22} Y F ‘L. V ’10. I '4 If" Tf . 1" { ml 7.13 3 15"“13 . "1.3436,". 3'; b3" ”13.1)3 '1I’3‘, ,IIA 1'. I gang. 1". .1 "'1 - 311"“ 1--.': .. 1914-2»; 'I3 1" Ii 2““ .14 ‘ .Url 3‘31"“ (‘9'... ....14'!‘:.'1"~3‘\‘ ‘i‘u F ”1.533 ' 3k 11:.1 1- w 1- 1 '1. - .1 2: ' {2. 3 ‘ I“ If. . ' '1 ‘3‘1 '33!IT31'!,1~':_.’ 33’ "31%). -I 3“ 3 i' 3 . 1.4.“ -‘ I '"I‘. l"l 3 V. 4 '. I “11" “ f': 3. .1 [1103?3'1'1“ 1h: {ltd-.1 :I'W “1““: ”"f‘ ' {'3 In O J‘H- v ..., v1. .10 , "I? 1 4-10 ">5'~1‘,II‘I “‘4’: 1y}: . ‘4‘”:‘S. 4. v (w '31-].1',"‘I'3 '1'} ‘; ’. .3, :. 't} 1.319% h._ : ‘3 ' ‘1‘“ ‘ “JP 335 1";IXLX' magi? r ..‘u 1 .‘1. ‘ l" 1 . ”I'I' I} ”h 33v'x'i ' 1153.. "‘{Zlg’ . ‘,?-‘ "-:'»-1 " 1'1 .3 . If.‘ - g 1:: » 41-" . I, ll,” 3’1“ | I . “3 s 'W ML 4.5 1;! i' 1 {1.111% ”9:: Ifldé‘c'ffi‘g 6'13th 4.» :"3'" ' .4 . "13??! . I 1 n.1,“: 1' . 0.2.1; I: f: H 1 . 17%.“ ~ ,3. #::3‘\:3, i‘fil :i?“~¢h§ In LA :1 , | ..I 1 {13'3“ 0m; ”a 3 ‘~1 . 1.: ‘4' 4 ...1 3111-" 11." .1 1w. 16.13.5131'1343' " 1793-1; {3: -.§‘13"rrr.~ 31‘1“ , 1 . .32, ., 1 - 1"‘(7'3'1‘3‘ 1.3:".‘ ...1'.“ I "M ' Mir-,2 1333-],' 5.713 ‘.‘ ' 4 ' 0 ”kn”? 3' 3 "'l‘ |1‘" '4 M511"?- ‘fl'f'; "1;: S I ‘I. I. .. "' ' "" 3 II" x'd" 31,3213. 1’ 1‘1 ",vafil' 3": J, 3' '1‘33,;J.‘%;[H}L;I:Cvm ' gas: I. khaki“ . ,. .‘MQO 1 , H '0‘ ,HI'I H 1 w i .1: PMI '4:'1 ' LT .- . ‘ Tun: “9.10": ‘7 JO 4 _\ D I‘D - I. . 41f ' In. '34 ‘is‘zr'l‘ I V‘I .,[.. a... . mfg-,1! , 1. . ...'.1.. 1,- I 1.1 , on. " L343 1.! ‘I"1‘.3j"1'31‘3 5.03.5 ,3-' «1-,. (L a" ‘ Kiy’2K‘H . .1 :91: ”£3.31‘33'9‘IA§3¥#|:.$ 1 \1 11 ‘5”; '1 '.|r"41{I.Hi' (‘3': $311 '15.?“ 8,3501 [51,211 P ._ ,J 1r ‘ "#1,“? 3" 1‘7 'r'I L .11‘611'41'11‘11fl-i; :10 :',.I‘A'3" In. 23.31‘J'Ul:¥3 i; g 1 “3391 333313-33“. 3 '-'. 11.411.391.513- , 1J7 .40" .1 ' :W x 1 I1' 5? '-~r'_ . .J_ 11" 1.. n :4: :13 7"” W ' 197'6131 |_ >r ‘ 1’3yj._,"$‘ 3E". 3?: "'_;1’ Skin" . ,WII‘L7 \3l ' 2'” v 3': all.“ I ‘1 . .4 . ‘1 , 31:1 " ' ‘ ... 1w: ‘. “I“ '1‘," ‘ .I' “Hutt ’1'!“ . 1'1.' 4‘: I1 ‘1‘! 1"!" . M: Tléh’ LI". '33. IIII.:'1|1'; 4:" » NIEV RS SIIYYLBRAR IIIIII III | gr: ~53, "'1'? dgi-wo ! IIIIIIII IIIIIIIIIIIIIIIIIIIIIIIIIIII prim. ram-2::- ‘zf' ’) “.~' :4 :3? ff; « a ‘. 7'3» 3 1293 00618 4943 I Tat-“... t. r :3: ‘ ' k’... My- .. ,. This is to certify that the thesis entitled ANANALYTICALANDEXPERIMENJRL INVESTIGATION OFTHE DYNAMCRESPCNSEOFAFOURBARMECHANISWCONSTRIDI‘ED MAWSCDELASTICCOMPOSITEDM'ERIAL presented by Wang Chun-hwa has been accepted towards fulfillment of the requirements for M. S. degree in Mechanical Eng. /4;\X 773A» Major professor 0.7639 MS U is an Affirmative Action/Equal Opportunity Institution MSU RETURNING MATERIALS: Place in book drop to remove this checkout from LIBRARIES .—_‘—. your record. FINES will be charged if book is returned after the date stamped below. . {)2 AN ANALYTICAL AND EXPERIMENTAL INVESTIGATION OF THE DYNAMIC RESPONSE OF A FOUR BAR MEaIANISM CONSTRUCI‘ED FROM A VISOOELASTIC COMPOSITE MATERIAL By Wang Chun-hwa Thesis Submitted to Michigan State University in partial fulfilment Of the requirement for the degree Of MASTER OF SCIENCE Department Of Mechanical Engineering 1984 ABSTRACT AN ANALYTICAL AND EXPERIMENTAL INVESTIGATON OF‘EHE DYNAMIC RESPONSE OF A FOUR BAR MECHANISM CONSTRUCTED FROM A.VISCOELASTIC COMPOSITE IATERIAL Wang Chun-hwa The intense competition in the international marketplace for mechanism systems which operate at higher speeds. are less noisy and more energy efficient than previous designs has imposed considerable preasures upon the machinery designer. This is because classical rigid-body analyses are unable to predict the elastodynamic phenomena associated with these new modes of operation. In order to respond to these commercial stimuli. mechanical systems need to incorporate members with high stiffness to weight ratios and also high strength to weight ratios. The work presented here develops appropriate finite element models for four- bar mechanisms constructed in elastic and viscoelastic materi- als. Experimental investigations into the effects of different link materials upon the dynamic flexural response of four bar mechanisms are also described. The correlations between the analytical and experimen- tal results for the midspan transverse deflections of the coupler and rocker links are good. thereby suggesting that these models may be used with confidence in the computer aided design of high-speed machine sys- tems fabricated in the commercial materials and also composite laminates. Acknowledgment The author wish to express his sincere thanks to Dr. B.S.Thonpson for his continual suggestions and guidance during the investigation. Sincere appreciation is extended to Dr. G.E.lase of the Department of letallurgy, Mechanics and Material Science and Dr. B.Fallahi of the Department of Mechanical Engineering for serving on the guidance commit- tee for their helpful suggestions. The efforts and assistance of Kr. C.K.Sung and other friends in Michigan State University are sincerely appreciated. Special credit must be given to every member of author's family for their encouragement and financial support that help to complete this work. TABLE OF CONTENTS List of figures V“ Chapter 1. Introduction 1 Chapter 2. Variational Theorem 5 2-1. Theoretical development 5 2-2. Finite element formulation 17 Chapter 3. Mathematical Model of a Flexible Linkage 28 3-1. Planar beam element 28 3-2. Shape functions of beam element 33 3-3. Element mass and stiffness matrices 36 3-4. Assembling the mass and stiffness matrices Of the mechanism system 39 3-5. Linkage model 41 3-6. Construction of the system matrices for the finite element model 42 3-7. Modeling the viscoelastic constitutive equations 44 Chapter 4. Experimental Investigation 49 4-1. Material characterization studies 49 4-2. Experimental apparatus 59 4-3. Instrumentation 51 Chapter 5. Experimental and Computational Results 65 Chapter 6. Discussion of Results 77 Chapter 7. Conclusions 30 81 References Appendex 92 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 2-1.1 2-2.1 2-2.2 3-1.1 3-2.1 3-4.1 3-5.1 3-6.1 3-7.1a 3-7.1b 3-7.2a 3-7.2b 3-7.3 3-7.4 LIST OF FIGURES Definition of Axis Systems and Position Vectors.---- 6 Deformations of the beam element 17 Finite difference representation— 22 A General Bemm Element. 29 The Deformations of a Beam Element. 34 The Relation between Global and Local Coordinate.-~--39 Four Bar Linkage Model with Flexible Coupler and Rocker. 42 Finite Element Model for the Simulation. 43 Spring. 46 Strain-Time Characteristic of Spring. 45 Dashpot. 45 Strain-Time Characteristic of Dashpot. 46 Kelvin Model. 48 Standard Linear Solid Model— 48 Result for Material Testing 51 Result for Creep Testing 52 Method to Obtain the Relaxation Function 54 Relxation Function Obtained from Material Testing--—-55 Relxation Function Obtained from Curve Fitting-----56 Natural Frequency of the Steel Link 58 Natural Frequency of the Composite Link 58 Joint connecting the flexible coupler and rocker----60 vii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure PAGE 2 4-9 Experimental Four Bar Mechanism 60 4-10 Experimental Appartus Flow Chart 52 5-1 Steel Coupler Response Operated at 342 RPM (Oscilloscope Photograph) 57 5-2 Steel Coupler Response Operated at 342 RPM (Computer plot for Both Results) 68 5-3 45 Composite Coupler Response operated at 280 RPM (Oscilloscope Photograph) 57 5-4 45 Composite Couper Response Operated at 280 RPM (Computer Plot for Both Results) 69 5-5 45 Composite Coupler Response at 205 RPM. -------- 71 5-6 45 Composite Rocker Response at 205 RPM. 72 5-7 45 Composite Coupler Response at 255 RPM. ---------- 73 5-8 45 Composite Rocker Response at 255 RPM. 74 5-9 45 Composite Coupler Response at 297 RPM. ---------- 75 5-10 45 Composite Rocker Response at 297 RPM. 76 Chapter 1 Introduction For the past decade. the design of machine members have always been based on the rigid body method of analysis. which means designers have traditionally considered all members of a mechanism to be rigid bodies [1-5]. With this rigid body approach. elastodynamic phenomena. such as dynamic stresses and vibrations. associated with link elasticity are neglected. This method can only be considered to be resonably accurate for those mechanisms which operate at low speeds. The major problem. from designer's point of view is that the traditional rigid-body ana- lysis and synthesis techniques are unable to predict the response of these mechanical systems with flexible members because mathematical models need to incorporate link flexibility. Recently. machines have been required to operate at higher speeds with more accurate performance characteristics. Because of the higher Operating speed creates increased inertial forces. the mechanism must be fabricated as a lightweight form-design to reduce the inertial loading. Unfortunately. mechanism with lighter members develop deformations and vibration due to external and internal forces; therefore. the perfor- mance of the mechanism might not be acceptable owing to these inaccuracies. In addition. failure of the mechanism might occuur if a vibration analysis has not been undertaken. Owing to this effect. the analysis of non-rigid mechanisms have became extremly important in the mechanism design area. Research in the field of mechanism design has progressed from studies of systems containing one [1-9]. or more. flexi- ble links [10-30] during the past few years as researchers have attempted to deveIOp viable mathematical models for designing high speed mechanism systems with high stiffness weight ratios and lighter weight members because classical techniques based upon rigid-body dynamics are unable to adequately predict the performance characteristics of these flexible systems. These kind of techniques would allow a designer to deveIOp a light weight. flexible mechanism that would meet both require- ments. Comprehensive review articles of the early research in the area of the dynmmic analysis of elastic mechanism have been presented in [31. 32]. In all of these references, the systems are assumed to be con- structed from homogeneous isotropic materials such as carbon steels or aluminum alloys. Some of the more recent work were devoted to deveIOp- ing Optimum lightweight form-designs for the menbers of linkages based on optimizing the stiffness characteristics or focused on the cross-sectional dimensions and shapes of the links [33-37]. Hence links were designed with tapers and complex cross-sectional shapes which directly increases the cost of manufacture. while reducing product marketability. An alternative philosophy with which to design a mechanism with high strength-to-weight ratio and stiffness-to-weight ratio components is to fabricate the links from a modern fiber-reinforced composite. As is well known, these materials have much higher strength and stiffness to weight ratios than the commercial materials such as steel or aluminum; furthermore. the composite also have high material damping. and good fatigue life. Although considerable fundamental research has been undertaken on determining the mechanical properties [38-43] and the response Of the composite materials treated as elastic members [44.45]. the literature is devoid of reports on the dynamic viscoelastic response of mechanism systems built using these materials. Therefore. the objective here is to establish guidelines for the design of linkage machinery in composite materials. A four bar mechanism was constructed by incorporating these different materials as the coupler and rocker links and treating them as the flexible parts. A finite element code was developed for the nonli- near elastic analysis Of flexible four bar linkages based on a variational formulation [46] and presented in [47]. A second code was developed for a four bar linkage fabricated in a £45 graphite epoxy lam- inate. The material characterization tests show that the £45 degree composite was a viscoelastic material. There have been several deriva- tions of variational theorems in linear viscoelasticity [48-56]. The experimental four bar mechanism which was fabricated using a :45 degrees composite is analyzed by first developing an appropriate variational principle for a general multi-body system fabricated from a linear viscoelastic material prior to formulating a model of constitutive equa- tions of the composite laminate and then solving the resulting equations of motion by the finite element method. This single functional expression with its associated varational equation provides a complete description of this class of mixed boundary value problem and is relevant to the design of all mechanisms. By per- mitting arbitrary independent variations of the stress. strain. displacement and velocity parameters. this equation yields the governing differential equations as well as the relevant boundary conditions. As an illustrative example.an approximate solution is sought for the response of the flexible four bar linkage by developing a displacement finite element model of the system. This mathematical model incoropo- rates one standard linear solid model to represent the material's constitutive equation. The equation of motion is solved by numerical integration. and the analytical and experimental results are presented. Chapter 2 Variational Theorem The objective of this chapter is to develop the equations governing the motion of mechanisms constructed in viscoelastic materials. The approach follows the developments of reference [57]. 2-1 Ihgoretiggl Dggglopgent The volume of a viscoelastic body is taken to be V and the total surface area S is the summation of two prescribed area Sd and so' The dynamic problems of this viscoelastic body V are considered. Let o-x-y-z be a set of Lagrangian coordinates fixed in the body in a refer- ence state with zero deformations. strains and stresses. and furthermore, it is also assumed that these parameters have been zero throughout the previous time t. Employing a Cartesian tensor notation. at time t a general point P in the continuum has the general position vector ri. which is defined as where roi: the component measured in the o-x-y-z frame with the position Figure 2—1.1 Definition of Axis Systems and Position Vectors vector of the origin Of the body axes relative to the origin Of the inertial frame. rri: the position vector of point P in the reference state relative to the origin of the body axes. “i: deformation displacement vector. in Figure 2-1.1. O-X-Y—Z is an inertial reference frame. The field equations for the linear theory of dynamic viscoelastici- ty for an body describing a general spatial motion relative to O—X-Y-Z are as follows: (i). pi=zoi+3i+eijkdj(rok+rtk+uk) (2'1.2) this equation is obtained by differentiate equation (2-1.l). and is the velocity rate of change of position statement; where pi : absolute velocity associated with ti °ijk: alternating tensor dj : angular velocity defining the rotation of the Lagrangian frame (~).: the time rate of change with respect to the moving coordinate frame (') : the absolute rate of change with respect to time. (ii). (iii). (iv). The boundary conditions for prescribed displacement and tractions are written as ui=ui on Sd Eigaijnjgzi on S (2-1.3) Sd: surface on which the prescribed displacements are imposed S : surface on which the prescribed tractions are imposed 0' gi: surface traction vector ( ): a prescribed quantity n-: the outward normal from the surface Sd Strain-displacement relation eij=(1/2)(ni’j+nj.i) (2-1.4) e..: Lagrangian strain tensor 1J Equation of equilibrium 0 : body force per unit volume Lagrangian stress tensor : mass density of the material : rate-of change of absolute valocity (v). Strainrstress relation in relaxation form oij=Gijklsk1to>+I§cijkl(t-c)[aeklttiIdc]dz (2-1.6) where Gijkl: Relaxation function oij : stress tensor the first term in the equation above represents the response at t=0, while the second term of the equation is the rate of change during a time interval dt. Equation (2-1.6) can also be written as “ij‘Gijn' den (2—1.7) where ‘ is the convolution form Equation (2-1.7) is the constitutive equation of the material In addition. the following energies are defined 10 kinetic energy: T=(1/2)ppi‘8ijpj potential energy: U=(l/2)Gijk1‘dakl‘deij The Objective now is to establish a variational theorem which con- tains equation (2-1.l) to (2-1.7) in the first variation of the functional. When the first variation is set equal to zero, and indepen- dent variations are permitted then these governing equations are the stationary conditions of the functional. The task of determining the stationary conditions of this functional expression may be achieved using a the Lagrange multiplier approach to introduce the constraints into the functional and generalise Hamilton's principle [24]. The functional may be written P ( ) - 4" v11; .dIeiJ (1,2) (ui'j'I’nj ' 1)] P (3) -~ -~ - + Vii .dtpi roi ni Oijk‘i.(rok+rrk+uk)]dv P — + salg')‘d(ui-ui)dsa]dt (2—1.s) The first variation may be generated using the standard rules of the variational calculus and this procedure involves utilizing the 11 divergence theorem and also the symmetric properties of the tensors. On establishing the first variation. this is then set equal to zero. yield- ing sJ=O =1 {:{jvdtsT-suldv +jsdzi‘d(8roi+8rri+bui)dsd "' Pvflgyfileij-(l/Z) (“Lfnj ' i)]¢1v +jvgxgfl)edlpi-roi-Ei-sijkzej(rok+rrk+uk)]dv + .86”? ’ «iii-aims, +dpvxgg).d[58ij-(1/2) (ani.j+6uj ' i)]dv P (3) _ *‘sali ‘d(-6ui)dSa] dt (2 1.9) Considering the first term in the equation (2-l.9). the varition of the kinetic energy T may be written as Mater/apps}.i ‘(a/api)[(1/2)pp1‘511p1]5pi =ppi‘6pi (2—1.10) Similarly. for the potential energy density. 12 50=(dU/Bdeij)6deij =(3/3daij)[(1/2)Gijk1‘deij‘dekllbdaij 'Gijkl'd'k1‘5d‘ij (2—1.11) Consider the third term in equation (2-l.9). since 5“i,j=5j.i' this may be rewritten ( ) I411} 'dlbsij-(1/2)(8ui'j+6uj'1)]dv =vai;)'d[5‘ijldv’vaij)‘d5“1,jdV (2—1.12) The second term in the equation (2-1.12) can be integrated by parts and upon applying the divergence theorem this leads to (a) 8 (1) _ (1) Ivlij ‘dfini’j [8111 nj ‘dbuidS Ivlij 'j‘dbuidv - ‘1’ (1) - (.) ‘Isdlii njcd8u1d8d+Isakij njtdbuidsa Iv‘ij,j‘d5nidv (2-1.13) hence. from equation (2-1.12) and (2-1.13). the third term in equation (2-1.9) can be rewritten as Ivlgg)‘dlbeij-(1/2) (Oni'j'tbnj ’ 1)]dV g (1) _ (1) (1) _ vaij .dbeijdv Islij nj‘dbnids'tjvlij'j‘dbnidv (2 1.14) 13 the fourth term in equation (2-1.9) may be written as (a) - ~ - ~ - vai ‘dlbpi Oroi Oui eijkbdj‘(rok+rrk+uk)]dv ‘I4‘i3)‘45P15V'I41i’)‘dl5?oi+531+eijk631*(rok+rrk+uk)lav (2-1.15) the second term of equation (2-1.15). may be integrated by parts over time and rewritten as t1 (3) ~ ~ g Iti[jvki ¢d[5roi+5ui+eijksdje(tokfrrk+uk)1dv1dt ... t0 (3) ~ ~ - tljvli .d[5roi+5“i+°1jk5‘j‘(rok+trkfuk)]dv dt _ (2) ~ ~ t1 _ [Jvii tdIOroi+8ui+eikadj0(rok+rrk+uk)NV]to (2 1,16) the second term in the right hand side of the equation (2-1.16) is zero because varoations at the extremes of the time interval are not permitted. Thus (a) -~-~- vai ‘dISpi 6:01 Gui eijkbdj‘(rok+rrkfuk)]dv =ngg’)‘dapidv-Iglg"'d[8r°i+6ni+eijkbej‘(rok+rrk+uk)ldv (2-1.17) 14 therefore. the first variation of the functional I may be written as 81=O =I§:[Iv5pi‘d‘PPi+1§”)dV+I;6x§;’SdteiJ-(1/2)(ni.j+uj,1)]av +Ivalga).dlpi—EOi-Ei-eijkaj.(rokfrrkfnk)]dV +186“? )ed (Hi-uiMSa-PLdggij sufiLGij k1 ‘denldv +Ivd8ui'[1§33j+l§')]dV+8droi'[I;Ag')dV+JsaEidsd] +sddjslj;eijk(rok+rrk+uk)i§"dv+jsad5ni'(-x§"-a§;’nj)aso] -Jsd(;i-ui)‘dbgidsd+jsad6ui‘(ii-1:3)nj)d86]dt (2—1.13) By permitting independent arbitrary variations in the system parameters. the Lagrange multipliers may be determined and written *i')"PPi (1) - ‘13 ‘Gijki‘d‘k1’“ij 15 (a) - - A1 = 8i (2 1.19) So the final form for GI becomes =I::[Ivboij’d[eij-(1l2) (ui'j-O-uj , i)Iav 'vaapi.d[pi—;Oi';i-°ijkaj.(rok+rrk+uk)lav ‘Isasgi.d(Ei-ni)dso+jvdbni.(°ij.j-péi)dv +IsadOUi‘(Eroijnj)dSa+[jsagidSa-I;ppidV]‘dbroi +[Isa;i°ijk(tok+rrk+uk)dsd+vaii°ijk(rok+rrk+nk)].d5‘j +.[V(aij—Gijkl‘d8kl).dssijdv1dt (2—1.20) Independent arbitrary variations of the deformation displacement. stress. strain. absolute velocity. and the kinematic parameters defining the rigid body equation of motion enable equation (2-1.20) to yield the field equations and boundary conditions for this class of dynmmic viscoelastic problem because each integral must be independently equal to zero. The above variational principle represents a generalization to 16 the theory of viscoelasticity Of an elastodynamic variational theorem presented in reference [15]. The characteristic equations obtained from equation (2-1.20) actu- rately define the dynamic viscoelastic problem for the design of mechanisms fabricated in materials with linear viscoelastic prOperties. In order to obtain viable solutions. simplifying assumptions must be introduced and inparticular the constitutive equation must be modeled. In order to establish an industrially orientated solution methodology. the finite element method was selected as the vehicle for generating approximate solutions. A number of approachs establishing finite ele- ment models of viscoelastic media have appeared in the literature .such as references [58. 59. 60]. l7 2;; Finite Element Formulgtion The variational equation of motion may be'employed as the basis for a variety of finite element models depending on the geometrical shape of the body being analysed. the type of deformation theory assumed to be apprOpriate. the information sought from the analysis. the accuracy of the model for the constitutive equations. and the assumption of whether the material is isotropic or orthotropic. A finite element model is developed herein for determing the flexural response of the beam-shaped links of planar linkage mechanisms deforming in the plane of the mechan- ism. These linkages are assumed to be fabricated from graphite-epoxy laminates with orthotropic prOperties. The developments are based on publications [61, 62, 63]. The link deformations are governed by the Euler-Bernoulli beam theory. shown as Figure 2-2.l 93‘I'\ u... I’5 Figure 2-2.1 Nodal Degrees of Freedom Describing the Deformations Of The Beam Element 18 and the transverse deformation may be written as w=[N']{U] (2-2.1) where {U}={U,.U,,U,,U‘}T Assuming that the member deforms on the axial and flexural mades. then the axial deformations is u=u.-zw.x =[Nn]{Un]-z[N"][U] (2—2.2) where tuni={u,,u,i [N]: the row vector containing the shape functions 2 : the spatial variable prependicular to the beam section in the plane of the mechanism. . : the spatial dervative with respective to the axial spatial variable. This symbol was also employed to denote the absolute time dervative (see equation 2-1.2) but confusion should not exist because in this section it is confined to Operations on the shape functions. The absolute velocity field p(x.t) is related to the nodal values in the 19 discretized field by p=[N][P} (2-2.3) where {P}={Pn1p ssssssss spu‘} [N]=[1-(x/L).l-(x/L).1.x/L.x/L.1] In order that equation (2-1.20) be employed as the basis for a fin- ite element analysis. this tensor expression must be reformulated in linear algebra format vaaij.d[eij“(1/2)(ni,j+u 'i)]dv=jgd{86]T'[[e]-[N'][U]]dv (2-2.4) j Ivbpi‘dp[pi-roi-ui+eijk‘dj(r°k+rrk+uk)]dv =j;dtaP1T-pttpi-INRI{pRJ-tuiifiilav (2-2.5, Isa581‘d‘35‘uiIdso=IsadISsITsi{fit-{Nita}iasa (2-2.6, Igdbui*(aij.j-pPi)dv=Igd{60}T‘[NlT‘([bliel-OINJIP})dv <2-2.7) Isadeui‘(Ereijnj)dsa=ISad{501T‘[N]T‘({El-{g})dsa (2-2.3) (IsazidSa-jgpfiidv)‘d6r°i=d{erolT‘(IsaizldSa-I;p[N](Pldv) (2—2.9) [Isa‘.° ij k(rok+rrk+uk)dsd+jvp.P° ij k.(rok+rrk+nk) ] .OOOj 20 =dststjsd(;.z)as+j;(?xpiidvi (2-2.10) Iv(°ij’Gijk1‘dek1).dbeijdv=Iv{5‘xx}T‘d({“xxl'IG}.dfexx})dv (2-2.11) hence the equation (2-1.20) can be rewritten as OJ=O =j§:[j;tsaiTsatte1—rnvliniiav +jgaptspiTsripi-[NR]{Pal-[Nirtilav +ISGIdOgIT‘({U}-[N]{U})dsa +Igdtean‘[N]T([DJIcl-plNliP})dv +Jsad{OU}T‘[N}T({El-{gl)dsa +d{8ro}T'(Jsa{E}dsa- vptuimdv) ”1535 ‘IJsd(;xE)dS+Iv(rxpP)dv] +I;dtbs]T‘({exxl-{Gl‘dtexx})dv]dt (2-2.12) Focussing attention upon equation (2-2.7). this may be written as 21 [vatsuiTtIN1T(tni{oi-ptNi{Pi)av =IvdISUIT'[NlT([D]{exxl-plN]{Ul-p[NR]{PRJ)dv (2-2.13) by utilising the equation (2-2.5). The first term of equation (2-2.13) may be subjected to integration-by-parts over x to yield I;dteu1TstNiT(a/ax)taxxiav =[IAdtemT‘[NthanidAL-Ldtsm'rs[N' Men} dv (2-2 .14, where [N']=(3/dx)[N] As stated in [62] consider an one dimensional finite element of length L.and cross sectional area A represented by a cartesian reference frame. If ox is the axial coordinate. then the unit vector on the end faces has the direction cosine 1=11. m=0, n=0 and surface tractions are {g}=i{oxx] ' (2-2.15) thus. equation (2-2.8) leads to Isad{601T'[N]T({El-{g})ds =IAdISUIT‘IU]T({E}t{axx})dA::g (2—2.16) The first term in the right hand side of equation (2-2.14) may be written as 22 O'ijltl Gijltll .._._—~-——_——o———'—~“ """Tijltol Figure 2-2.2 Finite difference representation 23 IAOISUIT‘[NlTicxxldAlx=L7IAd{OUIT‘[NJTchxldA]x=o upon combining there terms with equation (2-2.l6). terms cancel to give IAd{OU}T‘[N]T{;]dA (2—2.17) which is part of the equations of equilibrium. The second term on the right-hand side of equation (2-2.14) may be written as -j;aisu1Ter'1tax,iav The the standard linear solid. presented in reference [63]. will be employed to model the viscoelastic constitutive behavior. hence Oij+floij=2u(éij+aeij) (2-2.18) Using Figure 5-2.2. which is a finite-difference representation of the materials' constitutive curve relating stress and strain. the left hand side of the equation (2-2.18) may be written as =6ij(t,)+(Aoij/2) (2-2.19) 24 and, in addition. dij=(1/At)[oij(t1)-oij(t,)] =Ao[ij ) lAt (2—2 .20) Thus. the L.H.S. of the equation (2-2.18) becomes éij+aaij=maij lAt was“ (t,+(a/2) )Aoij =[(1/At)+(B/2)]Aoij+fioij(to) (2-2.21) Similarly. applying the above procedure to the strain. the right hand side of equation (2-2.18) may be written as 2p[[(1/At)+(o/2)]Aeij+aeij(to)] and hence equation (2-2.18) becomes [(2/At)+B]Aoij=2p[[(2/At)+a]Aaij+20eij(t,)]-28oij(t.) (2-2.22) Considering a one-dimensional model. equation (2-2.22) can be written as [oxx]=B(A[exx+2o{sxx(t.)])-C{oxx(to)] (2-2.23) 25 where Ae(2/At)+a Bszp/[(2/At)+fll C=2B/[(2/At)+B] where a. B. a are constants. exx(t,) and oxx(t,) are the strain and stress at the previous time step of the numerical algorithm. Since {exx]=[N']{Ul (2—2.24) the constitutive equation may be written as -Ivd{OU}T‘[N']T[AB[N'J{U}+2“B[3xx(t°)l-C{"xx(t°)l]dv =-j;d{aUITs[AE[N'ITIN'JIU]+2aB[N'thexx(t,))-CIN'ltaxx(t,>1]dv (2—2.2s) by defining [K]=I;AE[N']T[N']dv [M]=I;[N]Tb[N]av tuni=j;[N1T.[NR1av and combining equation (2—2.25) with (2-2.17) and (2-2.13) the resulting equations of motion are d{SU}T‘[-[I]{U}-2aB[N']T{sxx(t.)}+C[N'lTioxx(t.)}-[M]{U}-[MR]{PRI+ISG[N]T [EldA]=O .26 This may be written in a simpler form as [x1tur+iuii01=tfvi,coi-tuni{PR}+ISGIN'JTIEJaA where {fviscol=C[N'JTchx(t,)}-2aB[N'leaxx(t,)) Therefore. the final form of the variational equation of motion becomes 5J=0 =I::[Jvd{6exxlr‘[{exxl-[N']{U}]dv +IgdtbPlT‘pIIPI-[NR]{Pal-[lefilldv +Jsad{bng*({U]-[N]{U])dsa +d{60}T‘[[K]{U}+[M]{Ul-{f vi‘col+[MR]{Pkl-IsalNth;}ds]dv +d{5ro]T‘(IS {EldSa-IgplNJIPldv) 0' 27 +d63j ‘ [Isa (Ex?) dS+Iv (rxpP) dv] +Ivd{sgxx}Tt({oxx}-[G}‘d{exxl)dv]dt (2-2-25) Chapter 3 Mathematical Model of a Flexible Linkage The objective of this chapter is to develop the finite element equations for a general planar elastic linkage which deforms principally in the axial and bending modes. In the course of this develOpment. the nocal displacements. accelerations and shape functions experssions are derived. The stiffness matrix and mass matris are also presented. The model for the viscoelastic material is in the last section. Some of the material presented in this chapter is based on reference [61]. 3-1 Plggar ngm Element A general beam element is shown on the Fig. 3—1.1 in twO reference frames. The global frame O—X-Y and the local frmme o-x-y. The x-axis of the local reference frame is parallel to the beam element axis. The dotted lines represent the rigid body position of the beam element and the solid lines show its elastically deformed configuration. The elas- tic deformation Of the beam element could be discribed by six nodal displacements. here we denoted as uz. u,. u,. u‘. u,. u‘. These dis- placements are located at the deformed positions of the node A and node B. From the figure shown. the following relationships may be esta- blished in the global coordinates. 28 29 w I Figure 3-1.1 A General Beam Element 30 XAo-XAfuIcosO~u,sinO YAo=YA+u1sin6+u,cos0 (3-1.1) 9A.=0+n, Differentiating equations above with respect to time. the velocity and the acceleration at node A can be expressed in the global coordinate IA.=1Afu;cosO-u10sinO-n,sinO—uaécose YA.=IA+u1sin9+u1OcosO+n3cosO-u30sin9 (3-1.2) éA'Bé'Irt'l’ and IA.=XAffi:cosO-ZuIOsinO-ulézcosO-uIOsinO-fi’sinO-2u,écosO-u,O’sine-u,OcosO YA.=YAffi1sin0+2n1écos6-u103sin0+u1Ocos0+fi,eosO-2u,0sinO-uzé’cosO-uaOsinO 9A.!“é'i'fi3 (3.103) The absolute accelerations in equation (3-1.3) in global coordi- nates can be expressed in the local frame system with the following relations 31 gAngo 3080+YAp sine yAé-XA.sin0+YA.cosO (3-1.4) 6Ag=°+a. combining equations (3-1.3) and (3—1.4). the resulting equations 2A.-2A+a,-.,é*—2a,e.,a yA.=“A+a,-n,é’-2a,é+u,a (3—1.5) 9Afi=§+33 Equation (3-1.5) is the absolute accelerations of node A Of the beam element in the local frame. Apply the same procedure to node B EB.=iB+fi‘-u‘O’-2u,O-u,9 ?B.=§B+fi,-u,Oz-2u‘é+u4§ (3-1.6) swam Where 2A. yA. O. in. ya. 5 describe the rigid body motion of the beam element and are all kinematic quantities. Defining [n.il‘ combining (3-1 ..I #8 d! H! II ‘4 E 3as E 5! r1 e: rs fl! rs fl! re d: rs d! r‘ .5), (3’1.6)e this may be written as: 32 .Ind {8:1}: (3-1.7) in this section -2u30 -u,9 ‘ +2510 +u1O +0 +0 -2£,6 -u,O ~2n‘0 +u‘0 +0 +0 1 {u.}={6r}+(ul+{fin}+{fic}+{6t} where E : absolute acceleration (3-1.7) (3-1.8) (3-1.9) 33 6 : rigid body acceleration : elastic acceleration a: #3 : normal acceleration : Coriolis acceleration 636! t: tangential acceleration [in]. {ac}. {at} are the elastic and rigid body coupling terms. If the rigid body velocity and accelerations are small compared with those of the elastic nodal deflection. then in the equation (3-1.9) the pro- duct terms in vectors {ficl' [ac]. [it] are also small compared to the corresponding terms in {fir}+{fi]. The analytical model presented here did not incorporated theese coupling terms. 3;;_Shgpe Fggctions of Be m Elggggt The most widely used finite element approximation for representing the continuous function being studied. is the polynomial. Consider the beam element present in Figure 3-2.1.the shape functions for two kinds of deformations must be formulated. First. for the axial deformation.u1 and u‘ the choice of this shape function is a linear polynomial in x. 34 Figure 3-2.1 The Deformations of a Beam Element "2 U5 Us + Us l— ; . I A, E "i f “4 u(x.t)'61(;)u1(t)+d‘(;)u4(t) =a§+b (3-2 .1 ) where 61(;) and 6‘(;) are shape functions. By virture of fact that u(x.t) must be such that u(0.t)=u,(t)=b u(L.t)=u‘(t)=La+b the functions 61(;) and d‘(;) must satisfy boundary conditions 61(0)=l . 61(L)=0 d,(0)=o . A.(L)=1 therefore 6,=(L-§)IL A‘JIL Secondly. for the bending displacement field. the shape functions will 35 be described as a cubic polynomial in x. Bence w(x.t)=d,(;)u,(t)+d,(;)u,(t)+d,(;)u,(t)+6‘(;)u‘(t) =01+u,;4a,;’+u4;’ I (3-2.2) w(0.t)=u,(t)=o1 w(L.t)=u‘(t)=a:+a,L+a,L3+a‘L’ " (0.t)‘u,(t)=a, w'(L.t)=u‘(t)=o,+2a,LrI-3O‘L3 hence a,=(1/L‘)(-3a,—2L¢,+3¢,-L¢.) a.=(2/L’)(t,+2L¢,-¢,- the shape functions are. therefore a,=1-(3§'/L’)+<2I’/L') d,=(3;‘/L’)-(2;'/L') a,=x-<2;’/L)+(;’/L’) (3-2.4) 4.»;‘IME'IL’i The shape function of the beam element may; therefore. be written 36 d,(§)=1-'x'/L ¢,('£)=3<(L—;)/L)‘-2( (L-;)/L) ’ a, (Than—BIL)“ d. (§)-§/L (3-2 .5) 4, (3:3 (inf-2511.) ' a, (;)=(L-;) GIL) ' The transverse displacement. w(;.t). may be written as: w('x'.t)=¢,<;) ...,” w, G) 4., (0+4, (E) .u, (the. (I) ...“) (3-2.6) and the axial displacement u(;.t)=61(;)ou1(t)+dg(;)-u.(t) (3-2.7) ummmmm Assume that the beam element has a uniform cross section. then the local mass and stiffness matrices are represpectively [i] and [E], than 37 , . 1/3 0 13/35 0 11L/210 1.2/105 Symetric [Ek- pAL (3-3.1) 1/6 0 0 1/3 0 9/70 13L/420 0 13/35 \0 -13L/420 -L2/140 o —11L/210 L2/105 . . \ EA/L O 12EI/L3 Symmetric O 6EI/L2 4EI/L [k]= -EA/L o O EA/L O -12EI/L3 -6EI/L2 O 12EI/L3 . o ear/1.2 ZEI/L o -6EI/L2 4EI/L/ (3-3.2) When deriving the stiffness properties of a beam finite element. using the small strain theory. it is assumed that the transverse dis- placements are independent of the axial displacements or forces. In reality. however. a compressive axial force would tend to increase any transverse displacement of the beam. thus effectively decreasing the 38 transverse stiffness of the beam. While a tensile axial force would have the Opposite effect. This dependence of the stiffness matrix upon the axial loading is called geometric stiffening. and could become important in mechanism analysis where large axial forces are known to occur and also when the beam is very slender . One approximate method of including this effect is to calculate a geometric stiffness matrix [k5] based on large strain theory that would represent the coupling between the axial and transverse displacements. The geometric stiffness matrix for a beam element. which has been derived [46] and [64] is ' O o o o 6/5 L/lO -6/5 LllO [k§]=(P/L) L/lO 2L’/1s -L/10 -L’/3o (3-3.3) o o o o -6/5 -L/1o 6/5 -L/10 L L/lO -L’/3o -L/10 2L’I15 where F : the axial force in the element This matrix [k3] represents the change in transverse stiffness due to an axial force in the element. Tb include this coupling effect in the equation of motion. the geomatric stiffness matrix is simply added 39 to the element stiffness matrix [R]. ummmmm Russians 9.! the Rachael; m The finite element method generally requires the use of a global frame in order to assemble the element into a model for particular com- ponent. Therefore. a transfer matrix is necessary so as to transfer the stiffness matrix and mass matrix from the local frmm to the global frame before the matrices of the system can be assembled to form a model for the complete mechanism system in which the links have differing orienta- tions relative to the global frame. Consider the beam element shown as Fig.3-4.1 Figure 3-4.1 The Relation between Global and Local Coordinate The relations between two different coordinate are: 40 u1=U1cosO+U2sinO u2=-UlsinO+U2cosO u3=U3 (3-4.l) u4=UgcosO+U5sinO u5=-UgsinO+U5cosO “6:06 Thus. a transfer matrix may be written as: cosO sinO 0 O O O -sinO cosO O O 0 O O 0 l 0 0 O [R]= (3-4.2) O O 0 case sine O O 0 0 -sinO cosO O \ O 0 0 O O l and the global matrix In]. is related to the local mass matrix [H] by the expression [nHRIT [E] [R] (3-4.3) similarly [:14an if] [a] (34.4) 41 3;; Linkgge Model For the finite element analysis. the linkage is regarded as an instantaneous sturcture at every position. Consequently. the stiffness and mass matrices are different at each mechanism position. Generalized coordinates representing deflections are assigned to every joint permit- ting the members to deflect in the horizontal or vertical directions. As illustrated above. the following steps should be satisfied in the computer program. 1) An idealization of linkage stucture is needed. This will require selection Of the type and the size of the finite element to generate the system mesh. 2) The system-oriented element mass and stiffness matrices are generated for each e1 ement . 3) These element mass and stiffness matrices are assembled systematically to develOp the mass and stiffness matrices of the total linkage system. 4) Determination of unknown model displacements of the problem involves solving a system of coupled ordinary differential equations. These equations are obtained by using the equilbrium condition at the nodes. 42 Figure 3-5.1 Four Bar Linkage Model with Flexible Coupler Flexible -—-.-—' ‘*~ Rigid The linkage model is shown in Figure 3-5.1. the crank element is consi- dered to be rigid. while each flexible link. coupler and rocker. is divided into six elements. The element at both ends of the link. consi- dered as joints made of aluminum. is treated differently because of the different material prOperties. In Figure 3-5.2. system-oriented gener- alized displacements are labeled to describe the stuctural deformation Of the linkage as well as to maintain compatibility between the elements and nodes. For instance. at the joint between the coupler and the rock- er. U14 and U15 represents the nodal translations and U36 and U16 describe the rotational deflection of the coupler and the rocker at that point. A rigid connection between two elements. will be simulated by only one rotational displacement. 3:6 anstggction of system mgtriges In Figure 3-6.1. appropriate displacements are labeled on each retaining compatibility at the nodes. The system mass and stiffness 43 i8 MR .5 c‘ 17 Is. ,a Anus Figure 3-6.1 Finite Element Model for the Simulation 44 matrices of ith element are [mi]= m} [iii In). T .— i=1. 2' 3.0.00.0...12 3-7, Modeling 1h; Visgoelggtig Chastitgtiyg Egggtions The model representing a viscoelastic media is constructed to simu- late the experimental behavior Of viscoelastic material and involve differential equations relating strain. stress. and time. An important characteristic of these materials. especially whrn subjected to dynamic loading. is that they exhibited a time and rate dependence that is com- pletely absent in the constitutive relations of elastic materials. Although these types of material have the capacity to respond instan- taneously they also exhibit a delayed response. Thus the materials have the combined capacity of an elastic material to store energy and the capacity of a viscous material to dissipate energy. Materials with this type of behavior are termed viscoelastic materials and they have been the subject of several texts [48. 65-68]. the tOpic for chapters in standard reference texts [69-71]. and the objective of numerous papers [72-76]. For dynamic viscoelastic problems it is necessary to construct the entire solution without relying on the static elasticity results since after all. dynamic situations involve wave propagation phenomena. strain rate effects and attenuation characteristics. In this regard. 45 finite element solutions are proposed herein based on approximate forms for the constitutive equations because an exact solution requires the time history of the material to be known for all time. which is of course impractcal. The basic elements commonly used in the model representation are a spring and a dashpot. A spring. shown as Figure 3-7.1a. represents an elastic solid. and exhibits instantantaneous elastic strain and elastic recovery as Figure 3-7.1b. The equation relating stress and strain for a spring in time domain is o=Es or e=o/E (3-7.1) A dashpot. as Figure 3-7.2a. represent a viscous element. and exhibits irreversible creep and permanent set as Figure 3-7.2b. The differential equation relating stress and strain for a dashpot in time domain is o=pe (3-7.2) One of the combinations of the two basic elements is the Kelvin model. as Figure 3-7.3. which consists of a spring and dashpot in paral- 46 o—wvw—w E Figure 3-7.1a Spring. £9. E 0 to it Figure 3-7.2b Strain-Time Relation of Dashpot. O {o ;t Figure 3-7.1b Strain-Time Relation of Spring. 47 lel. At all times the elogation a of the two elements is the same. and the total force 6 will split into a; (for the spring) and a, (for the dashpot) in whichever way to make a the same. When applied to this model 01=E1a a,=pé (3-7.3) and from these two relations o=ol+oa =E,e+pé (3—7.3) Various combinations of the basic elements are also possible. The standard linear model. as Figure 3-7.4. representing a first order linear differential equation of stress and strain. This is a three parameter model which consists of a spring in serious with a Kelvin unit. The differential equation relating stress and strain for this model are (p1/E1)o+[1+(E1/E3)]o=ulé+Ela (3-7.4) 48 l A1 .1 Figure 3-7.3 Kelvin Model. '52 [.1 Figure 3-7.4 Standard Linear Solid Model Chapter 4 Experimental Investigation 5-1, Mgterial Characterization Studigg The experimental results concerning the material properties were deveIOped by Mr.C.K.Sung and Mr.J.Cuccio under the instructions of Dr. B.S.Thompson in August 1983 in Michigan State University [78]. The main causeof link deformations in a flexible mechanism is either the bending. or the flexural. mode and the associated deflection field is governed by the flexural rigidity. which is the product of the Young's modulus (E) of the material and the second moment of the cross-section area of the link (I). The materials chosen for the experimental work were a low carbon steel and a graphite-epoxy laminate with a symmetrical ply layup of £45 degrees relative to the longitudial axis of the link. The modulus of elasticity of the steel link specimens were obtained from supplier's data sheets. and the specimens were not subjected to mechanical testing. However. the composite materials have a greater variability Of mechanical prOperties. hence the characteristics of the lmminates need to be examined and quantified carefully. In order to 49 50 determine whether it is elastic. viscoelastic. or plastic. The purpose for the testing was to investigate the mechanicl pro- perties of the 145 degrees graphite-epoxy laminates. At first. the specimens were subjected to dynamic testing. which required the specimen to gain a prescribed maximum load 0.258 MPa over a range Of time inter- vals. Thus Figure 4-1 presents the results of the tests performed on the :45 degrees laminate and the maximum stress level was reached in 0.5. l, 10. 100 and 1000 seconds during the four tests. The results presented in Figure 4-1 suggest that the behavior of the material is certainly dependent on the rate of application of the loading. Thus. the 145 degrees laminate is a viscoelastic material. and the response curve implies that the constitutive relationship of strain and stress is nonlinear. In order to verify the deductions made from these test data. anoth- er test was undertaken to study the creep response of the materials. The results are presented in the Figure 4-2. The creep data in Figure 4-2 verifies that the :45 degrees composite is truly a viscoelastic material. The following method was adopted in order to incorporate the data from the material characterization studies into the mathematical model. The objective was to determine the relaxation function relating stress and time. Firstly. the maximum strain was measured on the response 51 I a: II E I I .I. .m L . .III ,M m 83. . . :32 I ma 2 I. ... 9333. 1:0qu Mom :33— » o-~omn§-m< 3595: “52:24 he“ 3:353 Baum: m? H... 3.5; mm a I mm m... g m . .-LIII]I.(I.. l: 8— 2:3 . mm: .8: 8 LI .. , .5. 1 {DENI'QUJL} 52 mmwuuofl noouu new «annex «It oummvm I It! AOI‘I] [u l :T I ]] I! 't'...][ «M93.-— . ..Sa 2 _. .. 2.53 § \.. .53 é... .53 8a . 8mm” 5... .mx n.o.~. L. {<5st u.~omn\a-m< 8433: 35:75 Etuétgm 3.58 m: « ...H.] m . .s\\~\m .2fl\.2H...ND.w UlEdflfiflihiZS 53 curve extreme left. This corresponded to a stress of 0.258 MPa. at 0.5 second after the load was initially applied. Then using this magnitude of strain.a horizontal line was measured from the point of load initiation on each of the response curve in Fig- ure 4-1 and a vertical prependicular line consructed until if intersected each response curve on the increasing load (upper) portion of the curve. This permitted the stress to be obtained and assuming that the rate of application of load was constant (it was programed to be constant on the MTS testing machine) this operation permitted a stress-time graph to be plotted. This is presented in Figure 4-3. Then using a PRIME 750 curve fitting software "CURVFIT" by changing the values of the parameters A. B in the equation I(t)=(A-a)e‘“t+s (4-1) ' and the data presentedin Figure 4-4. which is a curve of the standard linear solid model. the curve presented in Figure 4-5 was generated. By plotting the constants E1, E3. p1 of the relaxation modulus. [79]. can be obtained from the following procedures G(t)=Ele't/"1+(E1/E2)(l-e't/"1) =iE1-(E1/E2)Je’“t+(E1/E2) (4-2) 80 E1=A E2=A/B 54 stress stress at time t A Strain constant strain—.- 3183" 4'3 Method to Obtain the Relxation Function 55 .2533. 3.33.: Iowa @2538 dozen—E noducuuom Vie 0.33m fixmmv a 1111 / 1111 r111 1111 1111 {DE-immmm 56 9:qu ergo icon .3333 .8323: nemuauuom are one»; .00». mi... .8. I. I. 8.. In . b p p _ k I/ H rl 1 j cat-immune) 57 [11:1 [C After the dynamic and creep testing of the laminate is completed. the response data provides the fundations of several investigations. The :45 degrees laminate has a complex response because the stress-strain relation is not linear. and the gradient depends on the strain rate. Obviously. a wide range of different Young's modulus could be determined from the experimental data depended on the assumptions that considered to be necessary. Because a link of any mechanism exper- iences rapidly fluctuating stress levels. the response curve on the extreme left of Figure 4-1 was selected to provide the basis for the Young's modulus because this records the highest strain rate of any spe- cimcn. The approach was to draw a tangent to the response curve out the lower end of the upper portion of the curve in the region recording. the initial response immediately following load application. The effective Young's modulus for this viscoelastic material was calculated to be 3.143x10‘lbf/in3 while a mean value of 2.834x106 lbf/ina was calculated for the line joining the points defining the maximum and minimum stress levels. The final objective of the material characterization studies was to determine the material damping of the material. Each link specimen was clamped at one end to develop a cantilever configuration prior to 58 al Frequency of the Steel Link Natur Figure 4-6 a1 Frequency of the Composite Link- Natur Figure 4-7 59 deflecting and releasing the other end and recording the transient vibration. Besides obtaining the natural frequency of the link from the trace on the oscilloscope screen .simple logarithmic decrement calcula- tions were undertaken to calculate the damping ratio (t). The results are presented in Figure 4-5 for steel and Figure 4-6 for composite lami- nate. Thus it is evident that Of the laminate has a much higher damping ratio than the steel . and this property can be utilised in mechanism design to eliminate undesirable vibrations. 5:2 EsnL___rimental Apparatus Specimens were prepared to form matched pairs in the two link materials. At the end of each link. two clearance holes were drilled. These accommodated socket screws which clamped each specimen to the bearing housing and permitted the experimental four bar linkage to be constructed. The flexible links were fixed to the aluminum bearing housing by two socket screw at either ends with a flat plate. The small plate which was shown in Figure 4-7 ensured a smooth load' transfer between the principle components of each link. The experimental four bar linkage presented in Fig.4-7 was located on a large cast iron table which was bolted to the ground and wall of the labaratory. It had a rigid crank.with a link length of 63.5 mm (2.5 inches). while the lengths of the two flexible links. the coupler and rocker. were changeable. 60 Figure 4-8 Joint connecting the flexible coupler and rocker Mechanism Figure 4-9 Experimental Four Bar 61 The coupler and rocker links were supported on matched pairs of ball bearing (FAG R3 DB R12) of 0.25 inches in bore. Every bearing was preloaded using a Dresser torque limiting screw driver with :1 in-lbf preloading. This procedure is to ensure that bearing clearance is elim- inated. The impact loading associated with bearing cleariances would cause the links to have large defledtions. Nevertheless. where the bearings were assemblied with large axial preloading. the deflections of the linkages will be decayed. Owing to these affects. the torque limit device must be employed to accurately preload the joints. A 0.75 hp variable speed D-C motor (Dayton 22846) powered the crank through a 0.625 inch diameter shaft supported by a cast iron pillow box bearings. A 4 inch diameter fly wheel was keyed to the shaft thereby providing a large inertia to ensure a constant crank frequency. when operating in unison with the motor's speed controller. 1;; Instrumentation The instrumentation flow chart in the experimental work is shown in the Figure 4-9. The rated speed of the electric motor was measured to three decimal places by a HP 5314A universal counter which was actived by a digital-magnetic pickup. model 58423. by Electro corporation. These devices allowed the Operator to adjust the speed controller of the motor in order to achieve the desired speed. The experimental results present the variation of link deflections Visual monitorino ' ' ' ' - °f HP 5314A Crank speed (rpm) universal counter 62 a 6‘ Electra-Corp. pickup 12 vein 5. D. C. LiNks strain gaged Mechanism ' 60 tooth' D.C. motor spur gear Air Pox pickup 12 volts D.C.. I speed controller I ,__I=== system zero crank-angle configuration DEC .. DPD 11/03 Dul ‘9 T mi crocomputel'. Hicromeasurements system 2100 iiheatstone bridge and amplifier 1 dynamic strains post-processed experimental resui ts V Figure 4-10 -—-..,, ...-~ -..-...-.. .. ... Experimental Appartus Flow Chart: 63 with crank angle. Strain gages were bonded to the midspans of the coupler and rocker links of the four bar mechanism and shielded cables were used to reduce the effect of electromagnetic fields which can introduce spurious noise into the signals. In fact. noise from electro- magnetic fields and other sources. which was superimposed upon the strain gagesignals. was considered to be a major signal conditioning problem in this experimental work. A low pass filter was builted to eliminate the high frequency noise from the signal and the filter had a variable cut-off frequency. In order to relate the strain gage Signal to the configuration of the experimental mechanism. another transducer arrangement was esta- blished. A zero velocity digital pick up. Airpax 14-0001. was located so as to sense the bolt head at the end of the crank when the four-bar-linkage was in the position of zero-degree crank angle. This mechanism configuration signal and the output from the gages were either fed to the oscilloscope with a C-5C camera attachment for photographically recording the response. or to a digital data acquisi- tion system (a DEC PDP-11/03 microcomputer with 5 Mb hard disk). The BNC cables from the experimental apparatus were connected to a input-output module. This device had 16 analog-digital channels. 4 digital-analog channels and two schmidt triggers. Using the code developed for digital data acqusition. the flextural response signal was recorded from the zero crank angle position through 360 degrees by fir- ing on of the schmidt triggers. 64 Chapter 5 Experimental and Computational Results The analytical and experimental results are obtained by using the finite element method and the equations were solved by the Newmark method [77]. The [K8] matrix desiribed in the chapter 3 was included in the mechanism model developed for the elastic material. Each flexible link was devided into six elements.and the element at both ends of the link considered as bearing housing made of aluminum that need to input the different material prOperties and dimensions into the simulations. These analytical results are then compared with experimentally obtained strain data in order to verify the correlation between the analytical and experimental results. The main assumptions made formulating the computer model for the simulations are 1) All bearings were considered frictionless and without clearance. 2) Out-of-plane motions were disregarded. 3) Only small elastic deformations from the rigid body equlibrium position were assumed. 4) Gravitational acceleration was considered to be smaller than the elastodynamic accelerations. 65 66 5) The crank speed was assumed constant. 6) The natural frquency and damping ratios were calculated from experimental work. In Figure 5-1. it shows the dynamic response of the coupler taken by C-5C camera from T-912 oscilloscope. The material used was a low carbon steel. Figure 5-2 presents a comparison of the dynamic responses of the analytical and experimental investigations. It shows good core- lation between the two. The dimensions of the link and other necessary data were listed below. Link lengths: Ground 16 inches Crank 2.25 inches Coupler 12 inches Rocker 12 inches RPM of the mechanism: 342 RPM Cross sectional area of the flexible links: Width 0.75 inch (in the plane perpendicular to the mechanism) Depth 0.055 inch (in the plane of mechanism) Young's Modulus 30 x 10‘ psi bug—...- ... . . ._ Figure 5-1 Steel Coupler Response Operated at 342 RPM (Oscilloscope Photograph) Figure 5‘3 45 Camposite Coupler Response operated at 280 RPM (Oscilloscope Photograph) 68 con 2N 2mm man an uncommom Hoamsoo Hoopm «In oueumm .mwm¢0mo_ Om m402< 22(10 r m Haemosw. Hmown 00 OmquOP-Oz 69 2mm cam on uncommom .HOHmsoo opwmomsoo common n: vim one»; Ammmmumav “WEE 538 can a: an. no a and: _ > > .x P . >..:.o- < a d << A < . \ I ..o >E 9 .mEoEtooxm\ _8_;_2> << .: .Hmpnoswhomxm Hmofibafi 00.0 40.. ..—-~ c:&1&.;1&1c:E*PiC>2:un 73 00m a 2mm mmm pm ammomwom Hmfimsoo opwmomsoo common w: Ammmmummv mqoz<.xz > I // Hmpnoswhomxm Hmoapaaae< /\ /\ As «0 vd C)EJ&.#JEJL>E+H+C>2:UJ Chapter 6 Discussion of Results Figure 5-2 shows that the corelation between the analytical results and the experimental response are extremely good. thus proving that the variational theorem has been successfully applied and the correct model has been used to formulate the elastic linkages. However. the results for the composite material presented in the previous chapter do not have such a good correlation between the analytical and the experimental responses. The difference in amplitude between the two results at 280 RPM for the coupler is about 0.09 mm. which is an extremely small deflection. as shown in Figure 5-4. All the rocker responses have larger differences between two curves if compared with the coupler. Plausible explanations for the difference in the two classes of waveform will now be discussed. Firstly. the alignment of the experimental mechanism introduces an error. Since several specimens with different lengths were employed. upon completing each set of tests it was necessary reconstruct the mechanism : therefore. it was necessary to change the position of the rocker ground joint. Whenever this joint was moved back or forth. the alignment had to be carefully reconfigured. Furthermore. the clearance holes on the specimens for the socket screws will affect the alignment. The alignment can only be improved by considerable care prior to under- 77 78 taking the experimental work. Secondly.the value of the second moment of area of model for the joints in the mathematical representation are only approximate values. The cross section area of the real joints are rather complex. while the cross-section area being used in the simulation are the values of rec- tangular cross sections. Also. the lengths of the joints used in the simulation are the lengths from the centerline of the bearing to the other end of the joint. these lengths are a little bit shorter than the true lengths of the joints. However. the model of the joint did include the total mass of the joints. Thirdly. the dynamic responses of the experimental work were often disturbed by unknown sources of electromagnetic noise and the results are not always repeatable. Thus. it was necessary to take several data sets for the same RPM and select the best curve after first eliminating the noise disturbance using the Fast Fourier Transform program FILTER developed in the laboratory. This was developed using FORTRAN and MACRO. Further discussion concerning the stability Of the response in the experimental work has been discussed in [51] but is outside the scope of this thesis. Forthly. the sampling rate of the digital data acquisition may not be high enough to adequately represent the signals being monitored. In the present program . the PDP-ll can pick up one datum in 0.000363467 second: therefore. the curve obtained from PDP-ll were often more than 79 one revolution. If the experimental result is to be superimposed on the analytical result. judgements. based on the oscilloscope photograph were necessary in order to decide the starting point as well as the ending point of the revolution. Misjudgement of this activity may be resulted in a phase shift. The technique of calibration is considered to be the most critical reason which resulted in the difference between two response curves. This is because that the results of the calibration govern the magnitude of the experimental results. To do the calibration. a load was applied at the midpoint Of the specimen which was configured as a simply sup- ported beam. then recorded the deflection read from the dial gage. At the same time. the response was recorded from the oscilloscope so as to measure the voltage (strain) developed by the gage. The problem here is that the dial gage. must be observed by the operator in addition to recording the oscilloscope deflection and also apply the load on the link. Thus there are tOO many operations to be performed simutaneously. Two experimentalists. operating in unison improved the accuracy of this procedure. but a more accurate instrument for calibration in considered to be necessary. such as a calibration fixture. Chapter 7 Conclusions A variational theorem has been developed and has been shown to pro- vide a vaiable formulation for the finite element analysis of an experimental linkage fabricated with a viscoelastic material using a three parameter solid.model representing the material constitutive equa- tions. Although differences between the experimental and analytical results exist. these are rather small if attention is focused on the qusi-static response. The limitations in predicating the dynamic response are probably due to the first order model of viscoelastic material. As discussed previously. advanced approaches exist and these should probably be adopted a future work. The simulations do. however. provide a conservative prediction for the dynmnic response hence this model could be used with confidence in an industrial computer aided design environment for the design of high speed mechanism systems and robotic manipulators. 80 81 References 1. Smith. M.R. and Maunder. L.. "Inertia Force in a Four Bar Linkage." Journal of Mechanical Engineering Science. 1967. Vol. 9. NO. 3. pp.218-225. 2. Elliott. J.L. and Tesar. D.. "The Theory of TOrque Shaking Force. and Shaking Moment Balancing of Four Link Mechanism." ASME Jour- nal of Engineering for Industry. 1977. Vol. 99B. pp. 715-721. 3. Orlandea. N.. Chace. M.a. and Calahan. D.A.. "A Sparsity-Oriented Approcah to the Dynamic Analysis and Design of Mechan- ical Systems-Part 1 and 2.” ASME Journal of Engineering for Industry. Vol. 99. No. 3. 1977 . pp. 773-779. 4. Chace. M.A. and Smith. D.A.. "DAMN. A Digital Computer Program for the Dynamic analysis of Generalized Mechanical System." SAE Paper 710244. SAE Annual Meeting. Detroit. Jan. 1971. 5. Sheth. D.N. and Vicker. J.J.. "IMP(Intergrated Mechanism Pro- gram) A Cbmputer-Aided Design Analysis System for Mechanism and Linkages." ASME Journal of Engineering for Industry. Vo1.94. 1972. pp. 454-464. 82 6. Thompson. 8.8. and Barr. A.D.S.. "A Variational Principal for the Motion of Components of Elastic Mechanisms." 1975. Proceeding of the Fourth World Congress on the Theory of Machines and Mechanism. Newcastle Upon Tyne. U.K.. Paper No. 43. 7. Viscomi. B.V. and Ayre. R.S.. ”Nonlinear Dynamic Response of Elastic Slider-Crank Mechanism."ASME Journal of Engineering for Indus- try. 1971. Vol. 938. pp. 251-262. 8. Neubauer. A.H.. Cohen. R. and Hall. A.S.. ”An Analytical Study of the Dynamics of an Elastic Linkage." ASME Journal of Engineering for Industry. 1966. Vol. 88B. pp.311-317. 9. Jasinski. P.W.. Lee. H.C. and Sandor. G.N.. ”Vaberations of Elastic Connecting Rod of a High-Speed Slider-Crank Mechanism." ASME Journal of Engineering for Industry. 1971. Vol.94. pp. 636-644. 10. Sadler. J.P.. "On the Analytical lumped Mass Model of an Elas- tic Four-Bar Linkage Mechanism.” ASME Journal of Engineering for Industry. 1975. Vol. 97B. pp. 561-565. 11. Sadler. J.P. and Sander. G.N.. "A Lumped Parameter Approach to the Viberation and Stress Analysis of Elastic Linkages." ASME Journal of Engineering for Industry. 1973. Vol. 958. pp. 549-557. 12. Sadler. J.P. and Sander. G.N.. "Nonlinear Viberation Analysis 83 Of Elastic Four-Bar Linkages." ASME Journal of Engineering for Industry. 1974. Vol. 968. pp. 411-419. 13. Kalaycioglu. S. and Bagci. C.. "Determination of the Critical Operating Speeds of Planar Mechanisms by the Finite Element Method Using Planar Actual Line Elements and Lumped Mass Systems." ASME Journal of Mechanical Design. Vol. 101. NO. 2. 1979. pp. 210-223. 14. Bagci. C. and Kalaycioglu. 8.. "Elastodynamics of Planar Mechanisms Using Planar Actual Finite Line Elements. Lumped Mass Sys- tems. Matrix-Exponential Method. and the Method of Critical-Geometry-Kineto- Elasto-Statics." Asme Journal Of Mechanical Design. Vol. 101. No.3. 1979. pp. 417-427. 15. Thompson. 8.8. and Barr. A.D.S.. ” Varitional Principle for the Elastodynamic Motion of Planar Linkages." ASME Journal of Engineer- ing for Industry. VOl. 98B 1976. pp. 1306-1312. 16. Jandrasits. W.G. and Lowen. 6.0.. "The Elasto-Dynamic Behavior of a Counter-weighted Rocker Link with an Overhanging Endmass in a Four-bar Linkage. Part 1- Theory. Part II - Application and Experiment.” ASME Journal of Mechanical Design. Vol. 101. No. 1. 1979. pp. 77-98. 17. Imam. 1.. Sandor. G.N. and Kramer. S.N.. "Deflection and Stress analysis in High Speed Planar Mechanisms with Elastic Links." 84 ASME Journal of Engineering for Industry. 1973. Vol. 95B. pp. 541-548. 18. Erdman. A.G.. Sandor. G.N. and Oakberg. R.G.. "A General Method for Kineto—Elastodynamic Analysis of Mechanisms." ASME Journal of Engineering for Industry. 1972. Vol. 94B. pp. 1193-1205. 19. Midha. A.. Erdman. A.G. and Frohrib. D.A.. "An Approximate Method for the Dynamic Analysis of Elastic Linkages." ASME Journal of Engineering for Industry. 1977. Vol. 99B. pp. 449-455. 20. Midha. A.. Erdman. A.G. and Forhrib. D.A.. "A Closed-Form Numerical Algorith. for the Periodic Response of High-Speed Elastic Linkages.” ASME Journal of Mechanical Design. 1979. Vol. 101. No.1. pp. 154-162. 21. Midha. A.. Erdman. A.G. and Frohrib. D.A.. "A computationally Efficient Numerical Algorithm gor the Transient Response of High Speed Elastic Linkages.” ASME Journal of Mechanical Design. 1979. Vol. 101. No. 1. pp. 138-148. 22. Winfrey. R.C.. "Elastic Link Mechanisms Dynamics." ASME Jour- nal of Engineering for Industry. 1968. Vol. 91B. pp. 268-272. 23. Bahgat. B.M. and Willmert. K.D.. "Finite Element Viberational Analysis of Planar Mechanisms." Mechanism and Machine Theory. 1976. Vol. 11. pp. 47-71. 85 24. Thompson. 8.8.. "A Variational Formulation for the Finite Ele- ment Analysis of High-Speed Machinery." ASME Paper NO. 80-WA/DSC-38. 25. Allen.R.R.. "Dynamics of Mechanisms and Machine Systems in Accelerating Reference Frames." ASME Journal of Dynamic Systems. Meas- urement. and Control. Vol. 103. 1981. pp. 395-403. 26. Allen. R.R. and Dubowsky. S.. "Mechanisms as Cbmponents of Dynamic Systems: A Bond Graph Approach." ASME Journal of Engineering for Industry. Vol. 99. NO. 1. 1977. pp. 104-111. 27. Karnopp. D. and Margolis. D.. "Analysis and Simulation of Planar Mechanism Systems Using Bond Graphs." ASME Journal of Mechnical Design . Vol. 101. No. 2. 1979. pp. 187-191. 28. Tadjbakhsh. 1.0.. "Stability of Motion of Elastic Planar Link- ages with Application to Slider Crank Mechanism." ASME Journal of Mechanical Design. 1981. 29. Sanders. J.R. and Tesar. D.. " The Analytical and Esperimen- tal Evaluation of Viberatory Oscillations in Relaistically Proportional Mechanisms." ASME Journal of Mechanical Edsign. vol. 100. No. 3. 1978. pp. 762-768. 30. Yang. A.T.. Pennock. G.R. and Hsia. L.M.. ”Stress Fluctuation in High-Speed Mechanisms." ASME Journal of Mechanical Design. Vol. 103. 86 No.4. 1981. pp. 736-742. 31. Lowen. 0.6. and Jandrasits. W.G.. "Survey of Investigation into the Dynamic Behavior of Mechanisms Containing Links with Distrubut- ed Mass and Elasticity." Mechnism and Machine Theory. Vol. 7. 1972. pp.3-17. 32. Erdman. A.G. and Sandor. G.N.. "Kineto-Elasto-Dynmmics - A Review of the State-of-Art and Trends." Mechanism and Machine Theory. Vol. 7. 1972. pp. 19-33. 33. Imam. I. and Sandor. G.N.. "A General Method for Kinetoelas- todynamic Design of High-Speed Mechanisms." Mechanism and Machine Theory. Vol. 3. NO. 4. 1973. pp. 497-516. 34. Imam. I. and Sandor. G.N.. "High Speed Mechanism Design - A General Analytical Approach." Journal of Engineering for Industry. Trans. ASME. Vol. 97. 1975. pp. 609-628. 35. Khan. M.R.. Thornton. W.A. and Willmert. K.D.. "Optimality Criterion techniques Applied to Mechanism Design." Journal of Mechanical Design. Trans. ASME. Vol. 100. 1978. pp. 319-327. 36. Thornton. W.A.. Willmert. K.D. and Khan. M.R.. "Mechanism Optimization via Optimality Criterion Techniques." Journal of Mechanical Design. Trans. ASME. VOl. 1979. pp. 392-397. 87 37. Cleghorn. W.L.. Fenton. R.G. and Tabarrok. B.. "Optimimum Design of High Speed Flexible Mechanisms." Mechanism and Machine Theory. Vol. 16. NO. 4. 1981. pp. 399-406. 38. Jones. R.M..Mgnhnn19§ gt Cnmnostte Mntgrigls. McGraw-Hill book Company . 1 975 . 39. Tsai. S.W.. Halpin. J.C. and Pagano. N.J. (eds.). Egnposite Matgrialt Workghon. Technomac Publishing Company. Sanford. CA. 1968. 40. Christensen, R.M.. Mggnanigg gt gpnnntttg Mntetiglt. Wiley. NY. 1979 41. Sendeckyj. G.P. (Ed.). Mgghnnics gt Cgmnntite Mntgrintg. Com- posite Materials Series. Vol. 2. by Broutman. L.J. and Krock. R.M.. (Eds.) Academic. 1974. 42. Dietz. A.G.H. (Ed.). Cnmnositg Engingering Lnninntes. Wiley. NY. 1969. ‘3 a Borg. CaAa p charryg Fa]. ‘nd ElliOtt. SaYa (COOIdS a ) 9 Qomnositg Mntgrinls: Tetting nnn Degtgn. Thing gnntntnnng. March 1973. Williamsburg. VA. American Society for Testing and Materials. STP546. 44. Herrman. 0. (Ed.). Dynnmics ‘nt Sttugtured Saling. ASME Applied Mechanics Division Publication. 1972. 88 45. Lee. E.H. (Ed.). Dynamics _t Comnosite Materinls. ASME Applied Mechanics Division Publication. 1972. 46. Thompson. 8.8. and Sung. C.K.. ”A Variational Formulation for the Nonlinear Finite Element Analysis of Flexible Linkages. Part 1: Theortical DevelOpment." ASME Journal of Mechanisms. Transmissions and Automation in Design. submitted January 1984. 47. Sung. C.K.. Xing. T.M.. Wang. C.H.. Thompson. B.S. and Crow- ley. P.. " A Higher Order Variational Formulation for the Finite Element Analysis of Flexible Linkages. Part 11: Application and Experiment." ASME Journal of Mechanisms. Transmissions and Automation in Design, sub- mitted January 1984. 48. Christensen. R.M..Thgory gt Viscoelasticity; An Introdnction, Academic Press. NY. 1971. 49. Gurtin. M.E. and Sternberg. E.."0n the Linear Theory of Viscoelasticity.” Archives for Rational Mechanics and Analysis. Vol. 11. 1962. pp. 261-356. 50. Leitman. M.J.. "Variational Principles in the Linear Dynamic Theory of Viscoelasticity." Quarterly of Applied Mathematics. VOl. 24. 1966. pp. 37-46. 51. Schapery. R.A.. "0n the Time Dependence of Viscoelastic Varia- 89 tional Solutions." Quarterly of Applied Mathematics. Vol. 22. 1964. pp. 207-215. 52. Gurtin. M.E.."Variational Principles in the Linear Theory of Viscoelasticity." Archives for Rational Mechanics and Analysis. VOl. 13. 1963. pp. 179-191. 53. Leitman. N.J. and Fisher. G.M.C..”The Linear Theory of Viscoelasticity.” in Handvuch der Physik (S. Flugge. Ed.). Vol. 6a/3. Springer Verlag. Berlin. 1973. pp. 1-123. 54. Kline. R.A.."Variational Principles in the Linear Theory of Viscoelastic Media with Microstructure." Quarterly of Applied Mathemat- ics. Vo1. 28. NO. 1. 1970. pp. 69-80. 55.0nat. E.T.."On a Variational Principle in Linear Viscoelasticity." Journal of Mdcanique. Vol. I. NO. 2. June 1962. pp. 135-140. 56. Christensen. R.M.."Variational and Minimum Theorems for the Linear Theory of Viscoelasticity." Zeitschnift fur Angewandt Mathematik und Physik. Vol. 19. 1968. pp. 233-241. 57. Thompson. B.S. and Sung. C.K.. ”A Variational Theorem for the Viscoelastodynamic Analysis of High-Speed Machinery Fabricated from Com- posite Materials." ASME Paper No. 83-DET-6. 58. Zienkiewicz. O.C. and Watson. M.."Some Creep Effects in Stress Analysis with Particular Reference to Concrete Pressure Vessels.” Neuclear Engineering and Design. North-Holland Publishing Co.. 1966. pp 406-442. 59. Adey. R.A. and Brebbia. C.A.."Effecient Method for Solution of Viscoelastic Problems." Journal of The Engineering Mechanics Divi- sion. 1973. 60. Zienkiewicz. O.C.. Watson. M. and Kino. I.P.."A Numerical Method of Viscoelastic Stress Analysis." Interational Journal of Mechan- ical Science. Vol. 10. 1968. pp 807-827. 61. Midiha. A.. Erdman. A.G. and Frohrib. D.A.."Finite Element Approach to Mechanical Modeling of High-Speed Elastic Linkages." mechan- ical and Machine Theory. Vol. 13. 1978. pp 603-618. 62. Gandhi. M.V. and Thompson. B.S.. "The Finite Element Analysis of Flexible Components of Mechanical Sustems Using a Mixed Vatiational Principle.” ASME Paper NO. 80-DET-64. 63. Aoki. 8.. Iishimoto. K. Izumihara. Y. and Sakata. M..”Dynamic analysis of cracked linear viscoelastic solid by finite element method using singular element.” International Journal of Fracture. Vol.16. No.2 April. 1980. 91 64. Przemieniecki. J.S.. 1322;: 21 Mhiiil §I£l££l£11 Analysis. McGraw- Hill. 1971. 65. Flugge. W.A!1;ggnlnttjgttyt, Blaixdell Publishing Company. Waltham. MA. 1976. 66- Gross. Bu Minorities; Structure 9.: the Theories 91 thnnnlnntjnjtxt Hermann. Paris. 1953. 67. Bland. D.R-. Items): Linear Waiter. Perramon Press. Oxford. 1960. 68. Lockett. F.J.. Ngnltnnnt Viggoetggtig 521111. Academic Press. London. 1972. 69. Truesdell. C. and Noll. W.. in ‘gnnnhnnh fin; Phynit (S. Flugge. Ed.). Vol. 3. NO. 3. Springer Verlag. Berlin. 1965. 70- Mulvern. L-Eu Introduction Le the Mechanics 91. 11 222112212 Mgdjnm. Prentice-Hall. Englewood Cliffs. 1969. 72. Hashin. 2.. ”Viscoelastic Fiber Reinforced Materials." AIAA Journal. Vol. 4. No. 8. 1966. pp. 1411-1417. < 4‘06‘ 4'42 4.1‘ 73. Grot. R.A. and Achenbach. J.D.. "Linear Anisothermal Theory for a Viscoelastic Laminated Composite." Acta Mechanica. Vol. 9. 1970. Appendix 92 C00..0000.00.0.....000......0.00.00000.000.000.0000000000 CDCICDCDCDCDCSCSC5(DC)(5?)?!(Sflfl(QCSCDCDCDCDCSCDCDCSC) V III SSSSSSS CCCCCCC 0000000 V I S C C 0 O V I SSSSSSS C O O V V I S C C 0 O V III SSSSSSS CCCCCCC 0000000 THIS PRCXERAM IS DESIGNED TO SEEK m2 RESPONSE OF THE FOUR BAR MECHANISM CONSTRUCTED BY VISCOELASTIC COMPOSITE MATERIAL The mechanism have two flexible parts. coupler and rocker. Each flexible link consist of 6 elements. the two element at both ends are the joint elements that have different material properties and dimensions as the flexible link. The flexible part equaly devided into four element. We are seeking the step-by step solution by using the Newmark intergration method. (Ref: Bathe.K.J.and Wilson. E.L.." The Numerical Method in Finite Element Analysis" GGl.GGZ.NTA EL1.EL2.Gll(36.36) KLI(6.6).Q(36.1).MK(36.36).XG(40.40).KLJ(6.6) Gl(36.36).MLI(6.6).KX(40.40).MLJ(6.6) 00.0.00...0.00.00.00.00..........0.0.0.00000000000000000 ....OOGOQI’O‘IOOGDQ .0..000...0.0.0.000000000000000.0.0.0...0.00.00.00.00... YY(36.1).CC(36.36).OUT(36).'K(37.37).WORK(6.6) RA(36.1) .MASS .P(36.1) U(36).Ul(36).UZ(36).RCAP(36.36) Y(36.1).B(36.1).F(12) UL(36).OCAP(36).BB(36.1).NR(36.36) AC(40).KI(6.6).KI(6.6).KN(36.36) KL1(6.6).KL2(6.6).KL3(6.6).KL4(6.6) ML1(6.6).ML2(6.6).ML4(6.6).IL3(6.6) R4(6.6).RTH(6.6).R3(6.6).RT3(6.5) RI.RIl .RIZ.RIS .RI4.XLm1 .XLEVZJLDB .XLEM MASSl.MASSZ.IASS3.MASS4 XKLI(6.6).XKLJ(6.6).XKL1(6.6) XKL2(6.6).XKL3(6.6).XKL4(6.6) XMLI(6.6).XMLJ(6.6).XML1(6.6) XML2(6.6).XML3(6.6).XML4(6.6) BBl(1.6).BBZ(1.6).BBT1(6.1).BBT2(6.1) C cg C C THE DIMENSIONS FOR FINDING THE STRAIN AND STRESS O?O ('3 ll 93 ALPHA. BETA. DELT. MU . ETA EBSZ(1.1).EBS3(1.1).EBS4(1.1).IBSS(1.1) EBSS(1.1).EBS9(1.1).EBSIO(1.1).E3811(1.1) TAM2(1.1).TAM3(1.1).TA!4(1.1).TMM5(1.1) TAI8(1.1).TAI9(1.1).TAI10(1.1).TAI11(1.1) ARE mE STRAIN ATRESS AS mE INITIAL FOR NEXT STEP EBSBZ(1.1).EBSB3(1.1).IBSB4(1.1).EBSB5(1.1) EBSBS(1.1).IBSB9(1.1).EBSBlO(1.1).ERSBll(1.1) TAIB2(1.1).TAMB3(1.1).TAIB4(1.1).TA!BS(1.1) TllB8(1.1).TAIB9(1.1).TA!B10(1.1).TAMB11(1.1) U2L(6.1).U3L(6.1).U4L(6.1).USL(6.1) USL(6.1).U9L(6.1).UlOL(6.1).UllL(6.1) X2(6.1).X3(6.1).X4(6.1).X5(6.1) X8(6.1).X9(6.1).X10(6.1).X11(6.1) XT2(6.1).XT3(6.1).XTM(6.1).XT5(6.1) XT8(6.1).XT9(6.1).XT10(6.1).XT11(6.1) GS(36.1) Y2(6.1).Y3(6.1).Y4(6.1).Y5(6.1) Y8(6.1).Y9(6.1).Y10(6.1).Y11(6.1) YT2(6.1).YT3(6.1).YTA(6.1).YT5(6.1) YT8(6.1).YT9(6.1).YT10(6.1).YT11(6.1) O2(36.1) ZT2(6.1).ZT3(6.1).ZT1(6.1).ZT5(6.1) ZTB(6.1).ZT9(6.1).ZT10(6.1).ZT11(6.1) 0PEN( XL X6 61 ML MR ELl EL2 7.FILE='COM1') LOCAL STIFFNESS MATRIX GLOBAL STIFFNESS MATRIXAND MASS MATRIX WITHOUT B. C. GLOBAL STIFFNESS MATRIX INCORPORATING B.C.'S (1 TO 4) LOCAL MASS MATRIX GLOBAL MASS MATRIX INCORPORATING B.C.'S ELEMENT LENGTH FOR.THE COUPLER ELEMENT LENGTH FOR THE ROCKER XLEV (1 TO 4) LEVG'IHES FOR THE JOINTS E E8 YOUNG'S MODULUS OF EHE FLEXIBLE PART YOUNG'S MODULUS OF THE JOINTS RHO MASS DENSITY U Ul 02 DEFLECTION VELOCITY ACCELERATION OOWOGOOGOOOOOGOOOOO IN NEWMARX METHOD. 94 LOAD VECTOR ACC. VALUES FROM FLEX PRCBRAM CAP THE EFFECTIVE LOAD VECTOR BEFORE USING THE OED SUBRGJTINE 'LEOTlF' AFTERWARDS IT IS THE DYNAMIC DEFLECI‘IONS ME TO THE IDDE OF OPERATION OF THIS IMSL PACKAGE. THE DYNAMIC DEFLECI‘IONS IN THE LOCAL FRAME THE ACCELERATION VALUES FROM FLm AFTER THEY ARE CHANGED FROM RA. ”see INPUT 'IHE CONSTANTS IN THE CONSTITUTIVE mUATION OOVOOOOOGOOOOOOO PRINT 9. "IRE VALUE OF G1.G2.NTA' READ(1.‘)GG1.GGZ.NTA ALPHA=GGZINTA BETA=(GG1+GG2)/NTA MU=GGl/2. PRINT ‘. 'DJTER RPM OF THE MECHANISM' READ(1.‘)RPM “1?”? THE LDJG'IHES OF LINKS L1=12. L2=2.5 L3=9.0 L439.0 THE LENG'm OF EACH JOINT EEMENT XLE‘Jl=1 .5 XLENI2=1 .5 xuma=1 .75 XLEV4=1 .5 m1: LFNG'IHES OF THE M.EMEVT OF THE FLEXIBLE LINK GOO H.1=(L3-XLEV1-XLE‘12)/4 . m=(u-nm3-nm4 ) l4 . 9°W MATERIAL PROPERTIES FOR JOINTS M.E. mo=o.1 ES=10.3E+6 ?°? THE NECESSARY DATA FOR JOINT 1 ( ELE 1) WID1=O . 75 DEP1=O . 75 A1=WID19DEP1 ONE OF THE VISCO TERM (CONSIST OF THE STRESS TERM) ONE OF THE VISCO TERM (CONSIST OF THE STRAIN TERM) 95 MASSl=(RHO/384.)‘A1 RIl=(1.0/12.0)‘(WID1‘DEP1"3) (50!? THE NECESSARY DATA FOR JOINT 2 ( ELE 6) WID2=1.25 DEP2=0.75 A2=WID29DEP2 MASSZ=(RHO/384.)‘A2 RIZ=(1.0/12.0)‘(WID2‘DEP2“3) GO? THE NECESSARY DATA FOR JOINT 3 ( ELE 12) IID3=1.75 DEP3=O.27 A3=WID39DEP3 MASSS=(RHO/384.)‘A3 RI3=(1.O/12.O)‘(WID3‘DEP39‘3) THE NECESSARY DATA FOR JOINT’4 ( ELE 7) GOO WID4=1.25 DEP4=O.75 A4=WID49DEP4 MASS4=(RHO/384.)‘A4 RI4=(1.0/12.0)‘(WID4‘DEP4“3) an ('5 III MATERIAL PROPERTIES OF COMPOSITE(45 DEG) E:1.328E+06 DEP=O.135 RHOl=0.06 WID=1.36 AF'ID‘DEP MASS=(RHOl/384.) 'A RI=(l.O/12.0)‘(WID‘DEP"3) INITIALIZE THE CRANK ANGLE THZ=0.0 W“? 9°? INITIALIZE STRAIN AND STRESS FOR THE FIRST STEP IBSZ(1.1)=O. fll83(1.1)=0. EBS4(1.1)=O. IBSS(1.1)=O. EBSB(1.1)=O. fllS9(1.1)=O. EBSIO(1.1)=O. EBSll(1.1)=O. TAW2(1.1)=0. TAM3(1.1)=O. TAN4(1.1)=O. TAM5(1.1)=O. TA!8(1.1)=0. TAM9(1.1)=O. TAI10(1.1)=O. TAll1(1.1)=O. OGOOO INITIALIZE THE DEFLECI‘ION OF EACH EMT IN LOCAL FRAME THE BETWEEN U _L STANDS FOR THE ELEMENT N U2L(6.1)=0. U3L(6.1)=O. U4L(6.1)=O. U5L(6.1)=0. U8L(6.1)=0. U9L(6.1)=0. U10L(6.1)=O. UllL(6.1)=0. DEFINE TIME STEP AND RPM OF CRANK LINK C—g TH21=RPM‘2.0‘3.1415927/60.0 TSTEP=1.0 STEP=1.0/TSTEP ...? TO EVALUATE m5 CONSTANTS FOR THE NEWARK METHOD PAR=O.5 AA=0.25 Th1./(RPM/60.0‘360.0‘STEP) SSO=(1./(AA'T“2)) SSl=PAR/(AA‘T) SSZ=1./(AA‘T) SSB=(1./(2.‘AA))-1. SS4=PARIAA-1. SSS=(T/2.)‘(PAR/AA~2.) SSG=T‘(1.-PAR) SS7=PAR*T TO INITIALIZE U.UI.U2 AND UL AND O R“? 877 889 D0 877 I=1.36 U(I)=0. Ul(I)=O. U2(I)=0. UL(I)=O. Q(I.1)=O. CONTINUE DO 889 I=1.36 DO 889 J=1.36 KN(I.J)=O. 97 c: C BUILD UP THE MATRIX WITH SECOND DERIVATIVE C 0F SHAPE FUNC. 6 CALL SSP(EL1.BBT1.BBl) CALL SSP(EL2.BBT2.BBZ) C C DEFINE LOCAL STIFFNESS _MASS MATRICES C CALL LOL(E.A.EL1.RI.MASS.KLI.MLI) CALL LOL(E.A.EL2.RI.MASS.KLJ.MLJ) CALL LOL(ES.A1.XLEN1.R11.MASSI.KL1.ML1) CALL LOL(ES.A2.XLEN2.R12.MASSZ.KL2.ML2) CALL LOL(ES.A3.XLEN3.RI3.MASS3.KL3.ML3) CALL LOL(ES.A4.XLEN4.RI4.MASS4.KL4.ML4) C cssssessssseessssssssesassesssssssssssssssessesessss C C THIS IS THE START OF THE MAIN LOOP C csssssesssssssssseseesssssessssesssseessecesseessess DO 765 K=1.36O C C FIND THE ACCELERATION OF THE LINK C CALL KIN(RA.TH2.TH3.IH4.TH21.STEP.XLEN1. ‘ XLEN2.AC.XLEN3.XLEN4.L1.L2.L3.L4) C— C THE RIGHT'HAND SIDE OF THE EO. OF MOTION C IS -MK‘P(I). SO THE FOLLOWING STEP IS NECESSARY C; DO 990 I=1.36 990 P(I.1)=-RA(I.1) TT3=TH3 TTMBTH4 S3=SIN(TT3) C3-COS(T13) S4=SIN(TT4) C4=COS(TT4) BUILD UP TRANSFER MATRICES on? CALL RMTR(R3.RT3.TT3) CALL RMTR(R4.RTI.TTW) TRANSFER [K] MATRICES TO GLOBAL FRAME qnn CALL MMLT(WORK.KLI.R3.6.6.6) CALL MMLT(XKLI.RT3.WORK.6.6.6) CALL MMLT(WORK.KL1.R3.6.6.6) CALL MMLT(XKL1.RT3.WORK.6.6.6) 98 C CALL MMLT(WORK.KL2.R3.6.6.6) CALL MMLT(XKL2.RT3.WORK.6.6.6) C CALL MMLT(WORK.KLJ.R4.6.6.6) CALL MMLT(XKLJ.RT4.WORK.6.6.6) C C CALL MMLT(WORK.KL3.RM.6.6.6) CALL MMLT(XKL3.RT4.WORX.6.6.6) C CALL MMLT(WORK.KL4.R4.6.6.6) CALL MMLT(XKL4.RT4.WORK.6.6.6) (: C CONSTRUCT THE GLOBAL STIFFNESS MATRIX C CALL GLOLIN(XKLI.XKLJ.XKL1.XKL2.XKL3.XKL4.XG.GI) C C TRANSFER [M] MATRICES TO GLOBAL FRAM 6 CALL MMLT(WORK.MLI.R3.6.6.6) CALL MMLT(XMLI.RTB.WORK.6.6.6) C CALL MMLT(WORK.ML1.R3.6.6.6) CALL MMLT(XML1.RT3.WORK.6.6.6) C CALL MMLT(WORK.ML2.R3.6.6.6) CALL MMLT(XML2.RT3.WORK.6.6.6) C CALL MMLT(WORX.MLJ.RA.6.6.6) CALL MMLT(XMLJ.RT4.WORK.6.6.6) C C CALL MMLT(WORK.ML3.R4.6.6.6) CALL MMLT(XML3.RT4.WORK.6.6.6) C CALL MMLT(WORK.ML4.RA.6.6.6) CALL MMLT(XML4.RTM.WORK.6.6.6) C C CONSTRUCT IRE GLOBAL MASS MATRIX C CALL GLOLIN(XMLI.XMLJ.XML1.XML2.XML3.XML4.XG.MX) C IF(K.EQ.1)GOTO 10007 C (: C GANGE THE NOTATION OF THE DEF'LECI‘ION FORM C SYSTEM TO ELEMENTS (3: DO 201 I=1.6 J=I+1 201 U2L(I.1)=UL(J) DO 2011 181.6 2011 2012 2013 2014 2015 2016 2017 99 J=I+4 U3L(I.1)=UL(J) DO 2012 I=1.6 J=I+7 U4L(I.1)=UL(J) DO 2013 I=1.6 J=I+10 U5L(I.1)=UL(J) DO 2014 I=1.6 J=I+20 U8L(I.1)-UL(J) DO 2015 I=1.6 J=I+23 U9L(I.1)=UL(J) DO 2016 I=1.6 J=I+26 U10L(I.1)=UL(J) D0 2017 I=1.6 J=I+29 U11L(I.1)=UL(J) CALGJLATE THE STRAIN OF EACH FLEXIBLE ELE. COCO? THIS IS FOR THE EEMENTS ON THE COUPLER CALL CALL CALL CALL MMLT(I82.BBl.U2L.1 .6.1) MMLT(I83 .BBl .U3L.1 .6.1) MMLT(EBS4.BB1.U4L.1 .6.1) MMLT(I85 .BBl.U5L.1 .6.1) (500 1318 CALL CALL CALL CALL IS FOR THE ELEMENT ON m5 ROCKER MMLT(I88.BB2.U8L.1 .6.1) MMLT(EBS9.BB2.U9L.1 .6.1) MMLT(EB810.BBZ .U10L.1 .6 .1) MMLT(I811.BBZ .U11L.1 .6 .1) OOVO CALL CALL CALL CALL CALL CALL CALL CALL CALCULATE STRESS OF EACH FLEXIBLE EEMENT STRESS(TAW2 .TAWB2 . ISB2 .I82 .ALmA. BETA. MU) STRESS(TAW3 .TAWB3 . ISB3 .I83 .ALPHA. BETA. MU) STRESS(TAW4 .TAWB4 . ISB4 . I84 .ALPHA. BETA. MU) STRESS(TAW5 .TAWBS .ISBS .ISS .ALHIA. BETA. III) STRESS(TAW8 .TAWBS . EBSB8 . I88 . ALMA. BETA. m) STRESS(TAW9 .TAWB9 .ISB9 .IS9 .ALPHA, BETA. MU) 8TRESS(TAW10.TAWB10.ESB10.810.ALmA. BETA.MU) STRESS(TAWll .TAWB11 . EBSBll . I811 .ALPHA. BETA. MU) C C CALCULATE THE LOAD VECTOR Q=(M)"(P) 6 10007 c= CALL MMLT(Q.MX.P.36.36.1) 100 DELTb60./RPM/360. ETA=29DELT/(2+BETA‘DELT) C! C CALCULATE (B)‘TAW‘A ‘ETA DO 202 I=1.6 X2(I.1)=BBT1(1.1)‘T1l2(1.1)’A‘ETA X3(I.1)=BBT1(I.1)’TAW3(1.1)‘A‘ETA X4(I.1)=BBT1(I.1)‘TAW4(1.1)‘A‘ETA X5(I.1)=BBT1(I.1)‘TAW5(1.1)‘A'ETA X8(I.1)=BBT2(I.1)‘TAW8(1.1)‘A‘ETA X9(I.1)=BBT2(I.1)‘TAW9(1.1)‘A‘ETA X10(I.1)=BBT2(I.1)‘TAW10(1.1)‘A‘ETA X11(I.1)=BBT2(I.1)‘TAI11(1.1)‘A‘ETA CONTINUE TRANSFER X2....X5. X8....X11 TO GLOBAL FRAME DOC?O CALL CALL CALL CALL CALL CALL MMLT(XT2.RT3.X2.6.6.1) MMLT(XT3.RT3.X3.6.6.1) MMLT(XT4.RT3.X4.6.6.1) MMLT(XT5.RT3.X5.6.6.1) MMLT(XT8.RTI.X8.6.6 MMLT(XT9.RT4.X9.6. CALL CALL .1) 6.1) MMLT(XT10.RT4.X10.6.6. MMLT(XT11.RT4.X11.6.6. hfihl vv CALCULATE ONE OF THE VISCO TERM CALL GL05(XT2.XT3.XT4.XT5.XT8.XT9.XT10.XT11.GS) GOO CALCULATE (K)‘(U2L. . . . .UllL) CALL CALL CALL CALL CALL CALL CALL CALL MMLT(Y2.KLI.UZL.6 MMLT(Y3.KLI.U3L.6 MMLT(Y4.KLI.U4L.6 MMLT(Y5.KLI.USL.6 MMLT(Y8.KLJ.U8L.6 MMLT(Y9.KLJ.U9L.6 MMLT(Y10.KLJ.010L. MMLT(Y11.KLJ.UllL. GOO CALL CALL CALL CALL CALL CALL CALL CALL MMLT(YT2.RT3.Y2.6. MMLT(YT3.RT3.Y3.6. MMLT(YT4.RT3.Y4.5. MMLT(YT5.RT3.Y5.6. MMLT(YT8.RT4.Y8.6. MMLT(YT9.RTI.Y9.6. MMLT(YT10.RT4.Y10. MMLT(YT11.RT4.Y11. hihl vv 101 TIME ALPHA AND ETA ZT2(I.1)=YT2(I.1)‘ALPHA‘ETA/E 2T3(I.1)=YT3(I.1)‘ALPHA‘ETA/E ZTK(I.1)=YT4(I.1)‘ALPHA‘ETA/E ZT5(I.1)=YT5(1.1)‘ALPHA‘ETA/E ZT8(I.1)=YT8(I.1)'ALPHA'ETA/E ZT9(I.1)=YT9(I.1)‘ALPHA‘ETA/E ZT10(I.1)=YT10(I.1)‘ALPHA‘ETA/E ZT11(I.1)=YT11(I.1)‘ALPHA‘ETA/E ASSEIBLY ONE TERM OF THE VISCO PROPERTIES CALL GL06(ZT2.ZT3.ZT4.ZT5.ZT8.ZT9.ZT10.ZT11.0R) To CREATE TEE TERMS INSIDE TEE SECOND BRACKET ON LINE Bl P323 B! DO 890 I=1.36 Y(I.1)=880‘U(I) +882‘U1(I) +SS3‘U2(I) TO CREATE THE SECOND TERM ON LINE Bl M‘(AO‘UTWA2'U1T4A3‘U2T) CALL MMLT(B.MK.Y.36.36.1) TO CREATE LINE B1 P232 B! OCAP ETC THE EFFECTIVE LOAD QCAP(I)=Q(I.1)+B(I.1)+G§(I.1)-02(I.1) To CREATE mE EFFECTIVE STIFFNESS MATRIX LINE A4 9232 m1 RCAP(I.J)=Gll(I.J)+SSO‘MK(I.J)+881‘CC(I.J) TO SOLVE THE LINEAR EQUATIONS BY GAUSS ELIMINATION. CALL LINEO(OUT.OCAP.RCAP.WK.36.37.IERR) TO CREATE LINES B3 AND LINES B4 P232 B! G C LL DO 203 I=1.6 203 CONTINUE C: C C: (A C C C 890 C C C C C C C C DO 894 I=1.36 894 C C C LL DO 888 I=1.36 DO 888 I=1.36 888 C C (b C C ck DO 896 I=1.36 UD=U(I) UV=U1(I) UAFUZ(I) U(I)=OUT(I) 102 ’ U2(I )=SSO‘(U(I)-UD)-882‘UV-SS3 ‘UA U1(I)=UV+886‘UA+SS7‘U2 ( I) CONTINUE TREAT STRAIN AND STRESS AS WE INITIAL OF NEXT STEP one: ISBZ(1.1)=E82(1.1) TAWB2(1.1)=TAW2(1.1) ISB3(1.1)=IS3(1.1) TAWB3(1.1)-=TAW3(1.1) ISB4(1.1)=EBS4(1.1) TAWB4(1.1)=TAW4(1.1) ISBS(1.1)=I85(1.1) TAWBS (1.1)8TAW5(1.1) IISBS(1.1)=S8(1.1) TAWB8(1.1)-TAW8(1.1) IISB9(1.1)=89(1.1) TAWB9(1.1)=TAW9(1.1) EBSBlO(1.1)=I810(1.1) TAWB10(1.1)=TAW10(1.1) II$11(1.1)=811(1.1) TAWB11(1.1)=TAW11(1.1) CPOGG TO DEFINE mE DYNAMIC DEFLECTIONS IN THE LOCAL FRAMES 100 200 300 400 500 600 0 DO 100 I=2 .17 .3 J=I+1 UL( I)=(XJT(I)‘C3+(IJT(J) ‘83 DO 200 I=3 .18.3 J=I-1 UL(I)=-(XIT(J) ‘S3+(lJT(I)‘C3 DO 300 I=21 .33 .3 J=I+1 UL(I )=WT(I)‘C4+WT(J) ’84 DO 400 I=22 .34.3 J=I-1 UL(I )FWT(J)‘84+(IIT(I)‘C4 DO 500 I=1.19.3 UL(I)=(IIT(I) DO 600 I=20.35.3 UL(I)‘(IJT(I) UL(36)=WT(17)‘C4+(IJT(18)'S4 RK11=(K-1) /1 .0 0000 THIS IS FOR THE PLOT'TING . IT GIVES THE CRANK ANGLE. UL(9) DEFORMATION AT MIDPOINT 0F COUPLER 103 C UL(28) DEFORMATION AT MIDPOINT OF ROCKER WRITE(1.153)RK11.UL(9).UL(28) 765 CONTINUE 153 FORMAT(F14.8.1X.F14.10) CLOSE(7.8TATUS='KEEP') STOP END .iu.£u.s 2 104 COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO00.00....OOOOOOOOOOOOOOOOOtit L A 8888888 TTTTTTT OOOOOO L A A 8 T O O EEEEEE L A A 8888888 T O O L AAAAAAA S T O O EEEEEEE LLLLLLL A A 8888888 T OOOOOO OOOOOOOOOOOOOOOOOOOOOBO0.0....0......BOBO0.0000000000000000000000... This program is designed to seek the dynamic response of the elastic four bar mechanism. the material prOperties and the link length are changable. OOOOOOOOOOOOOOOOOO0......OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.00.0000. REAL EL1.EL2 REAL KLI(6.6).O(36.1).MK(36.36).XLJ(6.6) REAL Gl(36.36).MLI(6.6).MLJ(6.6) REAL YY(36.1).CC(36.36).OUT(36).WK(37.37).WORK(6.6) REAL RA(36.1) .MASS .P(36.1) REAL U(36).U1(36).U2(36).RCAP(36.36) REAL Y(36.1).B(36.1).F(12) REAL UL(36).OCAP(36).BB(36.1).MR(36.36) REAL KN(36.36).KG1(6.6).KG2(6.6).KG3(6.6).KG4(6.6).KI(6.6).KJ(6.6) REAL KGS(6.6).KG6(6.6).KG7(6.6).KG8(6.6).KG9(6.6).AC(40) REAL KGlO(6.6).KG11(6.6).KG12(6.6) REAL KGI(6.6).KGJ(6.6).KGM(6.6).KGK(6.6) REAL KS(36.36).KL1(6.6).KL2(6.6).KL3(6.6).KL4(6.6) REAL AST(2).ML1(6.6).ML2(6.6).ML4(6.6).ML3(6.6).R3(6.6).RT3(6.6) REAL RI.RIl.R12.RI3.RI4.XLEN1.XLEN2.XLEN3.XLEN4 REAL MASSI.MA882.MASS3 .MASS4 REAL XKLI(6.6).XKLJ(6.6).XKL1(6.6).XKL2(6.6).XKL3(6.6).XKL4(6.6) REAL XMLI(6.6).XMLJ(6.6).XML1(6.6).XML2(6.6).XML3(6.6).XML4(6.6) REAL KM1(6.6).KM2(6.6).KM3(6.6).KM4(6.6).KM5(6.6).KM6(6.6) REAL KM7(6.6).KM8(6.6).KM9(6.6).KM10(6.6).KM11(6.6).KM12(6.6) csssesss C C C C C C C C C C c c C C KL LOCAL STIFFNESS MATRIX XG GLOBAL STIFFNESS MATRIXAND MASS MATRIX.WITHOUT B. C. Gl GLOBAL STIFFNESS MATRIX INCORPORATING B.C.'S ML (1 TO 4) LOCAL MASS MATRIX MK GLOBAL MASS MATRIX INCORPORATTNG B.C.'S EL ELEMENTLENGTH FOR THE LINKAGE (COUPLER AND ROCKER) XLHV (1 TO 4) LEIGH FOR mE JOINT E YOUNGS MODULUS FOR.THE LINK DESIRED ET YOUNGS MODULUS FOR THE JOINTS RHO MASS DENSITY U DEFLECTION Ul VELOCITY U2 ACCELERATION (5(5t5t5t5t5t5t5(fitfitfitfitfi cessssss C IN NEWMARK METHOD. ceessess 105 C Q LOAD VECTOR C RA ACC. VALUES FROM FLEX PROGRAM C QCAP IS THE EFFECTIVE LOAD VECTOR BEFORE USING THE C SUBROUTINE 'LEQTIF' AFTERWARDS IT IS THE DYNAMIC C DEFLECTIONS DUE TO THE IDDE OF OPERATION OF C THIS IMSL PACKAGE. . C UL ARE THE DYNAMIC DEFLECTIONS IN THE LOCAL FRAME C P IS THE ACCELERATION VALUES FROM FLEX AFTER THEY C ARE CHANGED FROM RA. C: C INPUT THE LINK LmG'm DESIRI CT PRINT ‘.'WHAT KIND OF LINK LENGTH DO YOU WANT' PRINT ‘.'ENTER 1 FOR LONG LINKS' PRINT ’.'ENTER 2 FOR SHORT’LINKS' READ(1..)KLL IF(KLL.EQ.1)THEN GOTO 1 ELSE IF(XLL.EQ.2) THEN GO TO 2 ELSE GOTO 3 ENDIF C==s===== C THE LENGTHES OF THE LONG LINKS 1 L1=16.12 L2=2.25 L3=12.0 L4=12. C======= C THE LENGTHES OF THEISHORT LINKS 2 L1=12. L2=2.25 L3=9. L4=9. 3 CONTINUE C: C DECIDE WHAT KIND OF MATERIAL YOU WANT & PRINT ‘.'WHAT KIND OF MATERIAL YOU WANT!’ PRINT ‘.' 1. STEEL' PRINT ‘. ' 2. ALUMINUM' PRINT ‘.' 3. UNICOMPOSTE' READ (1.")XMAT IF(XMAT.m.1)mm GO TO 11 ELSEIF(XMAT. 190.2 )THEV GO TO 12 ESE IF(XMAT.EO.3)THEN GO TO 13 ESE GO TO 14 DID IF 106 DEFINE THE MATERIAL PROPERTIES MATERIAL PROPERTIES FOR ALUMINUM E:10.3E&06 WID=O.58 DEP=0.0795 A-WID'DEP RHO=0.100 DR-0.0048 XF-GO. BELTA=2.0‘DR‘XF OPEN(7.FILE='A') GO TO 14 C======= C 11 MATERIAL PROPERTIES FOR STEEL EF30.0EH06 WID=0.75 DEP=0.062 RHO=0.3 A=WID‘DEP DR=0.01 XF=31.1 BELTA§2.0‘DR‘XF OPEN(7.FILE:'ST') GO TO 14 C 13 MATERIAL PROPERTIES FOR UNIDIRECTIONAL COMPOSITE E320.5E+6 RHO=0.06 WID=O.75 DEP=0.08 AsWID'DEP DR=0.00984 XF=100.0 BELTA=2.0‘DR’0.0l/XF OPEN(7.FILE='UNIC') GO TO 14 CONTINUE n9: THE LENGTH OF EACH JOINT ELEMENT EJ=10.3E+6 RHOJ=0.1 XLEN1=1.5 XLEN2=1.5 XLEN3=1.75 XLEN4=1.5 0:? THE,LENGTH OF THE FLEXIBLE ELEMENT EL1=(L3-XLEN1-XLEN2)/4. EL2=(L4-XLEN3-XLEN4)/4. 0? THE DIMENSIONS FOR JOINT 1 ( ELE 1) 107 WID1=0.75 DEP1=0.75 A1=WID19DEP1 MA881=(RHOJ/384.)‘A1 RIl=(1.0/12.0)‘(WID1‘DEP1“3) THE DIMENSIONS FOR JOINT 2 ( ELE 6) WID2=1.25 DEP2=0.75 A2=WID2‘DEP2 MA882=(RHOJ/384.)‘A2 R12=(1.0/12.0)‘(WID2‘DEP2“3) on THE DIMENSIONS FOR JOINT 3 ( ELE 12) WID3=1.75 DEP3=0.27 A3=WID3‘DEP3 MASS3=(RHOJ/384.)‘A3 RI3=(1.0/12.0)‘(WID3‘DEP3“3) THE DIMENSIONS FOR JOINT 4 ( ELE 7) WID4=1.25 DEP4=0.75 A4=WID49DEP4 MASS4=(RHOJ/384.)‘A4 RI4=(1.0/12.0)‘(WID4‘DEP4‘*3) MASS=(RHO/384.) ‘A RI=(1.0/12.0)‘(WID‘DEP“3) INITIALIZE THE CRANK ANGLE TH2=0.0 DEFINE TIME STEP AND RPM OF CRANK LINK RPM=340.0 TH21=RPM‘2.0‘3.1415927/60.0 8TEP=1.0 TSTEP=1.0/STEP 03.? TO EVALUATE THE CONSTANTS FOR THE NEWARK METHOD PAR=0.5 AA=0.25 T51./(RPM/60.0‘360.0‘8TEP) 8808(1./(AA*T*‘2)) 881=PAR/(AA‘T) SS2=1./(AA‘T) 883=(1./(2.‘AA))-1. SS4=PARIAA-1. 885=(T/2.)‘(PAR/AA-2.) 886=T‘(1.-PAR) 887=PAR*T TO INITIALIZE U.U1.U2 AND UL AND O 108 D0 887 I=1.36 U(I) =0.0 U1(I)=0.0 U2(I)=0.0 UL(I)=0.0 O(I.1)=0.0 887 CONTINUE DO 889 I=1.36 DO 889 I=1.36 KN(I.J)=0. ”9: DEFINE LOCAL STIFFNESS _MASS MATRICES CALL LOL(E.A.EL1.RI.MASS.KLI.MLI) CALL LOL(E.A.EL2.RI.MASS.KLJ.MLJ) CALL LOL(EJ.A1.XLEN1.RIl.MA881.KL1.ML1) CALL LOL(EJ'.A2.XLEN2 .RI2.MASS2 .KL2.ML2) CALL LOL(EJ.A3.XLEN3.R13.MASS3.KL3.ML3) CALL LOL(EJ.A4.XLEN4.RI4.MASS4.KL4.ML4) cesssssessseessesssseesesssssesseessssssseseeesssss C THIS IS THE START OF'THEVMAIN LOOP ‘ cesesseesesssssassesseseeessssssssssseeesessssesess C CALCULATE THE ACCELERATION OF HE LINKS DO 765 K=1.360 CALL KIN1(RA.TH2.TH3.TH4.TH21.8TEP.XLEN1.XLEN2 * .AC.XLEN3.XLEN4.L1.L2.L3.L4) DO 990 I=1.36 990 P(I.1)=-RA(I.1) THIS IS NECESSARY BECAUSE THE R.R.S. OF THE EO. OF MOTION IS -MK‘P(I) GOO TT3=TH3 TT1=TH4 S3=SIN(TT3) C3=COS(TT3) 84=SIN(TT4) C4=COS(TTW) M? BUILD UP TRANSFER MATRICES CALL RMTR(R3.RT3.TT3) CALL RMTR(R4.RT4.TTH) M? TRANSFER [K] MATRICES TO GLOBAL FRAME CALL MMLT(WORK.KLI.R3.6.6.6) CALL MMLT(XKLI.RT3.WORK.6.6.6) CALL MMLT(WORK.KL1.R3.6.6.6) CALL MMLT(XKL1.RT3.WORK.6.6.6) CALL MMLT(WORK.KL2.R3.6.6.6) CALL MMLT(XKL2.RT3.WORK.6.6.6) CALL MMLT(WORK.KLJ.R4.6.6.6) CALL MMLT(XKLJ.RT4.WORK.6.6.6) CALL MMLT(WORK.KL3.R4.6.6.6) CALL MMLT(XKL3.RT4.WORK.6.6.6) CALL CALL CALL 109 MMLT(WORK.KL4.R4.6.6.6) MMLT(XKL4.RTK.WORK.6.6.6) GLOLIN(XKLI.XKLJ.XKL1.XKL2.XKL3.XKL4.XG.GI) TRANSFER [M] MATRICES TO GLOBAL FRAM CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL MMLT(WORK.MLI.R3.6.6.6) MMLT(XMLI.RT3.WORK.6.6.6) MMLT(WORK.ML1.R3.6.6.6) MMLT(XML1.RT3.WORK.6.6.6) MMLT(WORK.ML2.R3.6.6.6) MMLT(XML2.RT3.WORK.6.6.6) MMLT(WORK.MLJ.R4.6.6.6) MMLT(XMLJ.RT4.WORK.6.6.6) MMLT(WORK.ML3.R4.6.6.6) MMLT(XML3.RT4.WORK.6.6.6) MMLT(WORK.ML4.R4.6.6.6) MMLT(XML4.RT4.WORK.6.6.6) CALL GLOLIN(XMLI.XMLJ.XML1.XML2.XML3.XML4.XG.MK) IF(K.EO.1) GOTO 10007 DEFINE LOCAL NONLINER MATRICES FOR EACH KIND OF ELE CALL LONON(E.A.UR.EL1.UL.KI) CALL LONON(E.A.UR.EL2.UL.KJ) CALL LONON(EJ.A1.UR.XLEN1.UL.KGI) CALL LONON(EJ.A2.UR.XLEN2.UL.KGJ) CALL LONON(EJ.A3.UR.XLEN3.UL.KGK) CALL LONON(EJ.A4.UR.XLEN4.UL.EGM) CALCULATE KG FOR EVERY ELEMENT CON1=E‘A/EL1 F(1)=EJ‘A1/XLEN1‘UL(2) DO 101 I=2.5 DO 101 I=2.II.3 KK=J+3 F(I)=CON1‘(UL(KK)-UL(J)) CONTINUE F(6)=EJ'A2/XLEN2‘(UL(17)-UL(14)) F(7)=EJ‘A4/KLEN4‘UL(21) CON2=E‘A/EL2 DO 102 I=3.11 DO 102 J=24.33.3 M=J-3 F(I)=OON2‘(UL(J)-UL(M)) CONTINUE F(12)=EJ‘A3/KLEN3‘(UL(35)-UL(33)) DO 832 131 '6 DO 832 J=1.5 KG1(I.J)=F(1)IXLENI‘KGI(I.J) KG2(I.J)=F(2)/EL1‘KI(I.J) KG3(I.J)=F(3)/EL1‘KI(I.J) KG4(I.J)=F(4)/EL1‘KI(I.J) KG5(I.J)=F(5)/EL1‘KI(I.J) 101 102 110 KG6(I.J)=F(6)/XLEN2’KGJ(I.J) KG7(I.J)=F(7)/XLEN4‘KGM(I.J) KG8(I.J)=F(8)/EL2‘KJ(I.J) KG9(I.J)=F(9)/EL2‘KJ(I.J) KG10(I.J)=F(10)/EL2‘KJ(I.J) KG11(I.J)=F(11)IEL2‘KJ(I.J) KG12(I.J)=F(12)IXLEN3‘KGK(I.J) 832 CONTINUE C c: C C TRANSFER KG TO THE GLOBAL FRAME C CALL MMLT(WORK.KG1.R3.6.6.6) CALL MMLT(KM1.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG2.R3.6.6.6) CALL MMLT(KM2.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG3.R3.6.6.6) CALL MMLT(KM3.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG4.R3.6.6.6) CALL MMLT(KM4.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG5.R3.6.6.6) CALL MMLT(KM5.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG6.R3.6.6.6) CALL MMLT(KM6.RT3.WORK.6.6.6) C CALL MMLT(WORK.KG7.R4.6.6.6) CALL MMLT(KM7.RT4.WORK.6.6.6) C CALL MMLT(WORK.KG8.R4.6.6.6) CALL MMLT(KM8.RT4.WORK.6.6.6) C CALL MMLT(WORK.KG9.R4.6.6.6) C CALL MMLT(WORK.KG10.RM.6.6.6) CALL MMLT(KM10.RT4.WORK.6.6.6) C CALL MMLT(WORK.KG11.R4.6.6.6) CALL MMLT(KM11.RT4.WORK.6.6.6) C CALL MMLT(WORK.KG12.R4.6.6.6) CALL MMLT(KM12.RT4.WORK.6.6.6) 00090 CONSTRUCT THE GLOBAL [KGIMATRIX 111 CALL GLONON(KN.KM1.KM2.KM3.KM4.KM5.KM6.KM7.KM8. ‘ KM9.KM10.KM11.KM12.KX) TO CREATE mE LOAD VECTOR "Q". H00000 0007 CALL MMLT(Q.MK.P.36.36.1) TO CREATE THE TERMS INSIDE HE SECOND BRACKET ON LINE B1 P323 B! 000000 DO 890 I=1.36 890 Y(I.1)=SSO*U(I) +882‘U1(I) +883‘U2(I) TO CREATE THE SECOND TERM ON LINE B1 M‘(AO‘UT4A2‘U1T¥A3‘U2T) 000000 CALL M11189 MK: Ys36e36 '1) TO CREATE THE TERMS INSIDE THE.THIRD BRACKET ON LINE 00000 YY(I.1)=SS1‘U(I)+SS4‘U1(I)+SS§‘U2(I) 00 o O 0 DO 8901 I=1.36 DO 8901 I=1.36 CC(I.J)=BELTA‘MK(I.J) CALL MMLT(BB.CC.YY.36.36.1) cc ‘0 O H TO CREATE LINE B1 P232 B! OCAP ETC THE EFFECTIVE LOAD 000000 DO 894 I=1.36 OCAP(I)=O(I.1)+B(I.1)+BB(I.1) Q D & TO CREATE THE EFFECTIVE STIFFNESS MATRIX LINE A4 P232 B! 000000 DO 888 I=1.36 888 RCAP(I.J)=G1(I.J)+SSO‘MK(I.J)+KN(I.J)+SSI‘CC(I.J) 112 C C C TO SOLVE THE LINEAR EQUATIONS BY GAD SS ELIMINATION. C , CALL LINEQ(OUT.QCAP.RCAP.'K.36.38.IERR) C C: C C TO CREATE LINES B3 AND LINES 34 P232 H! C DO 895 I=1.36 UD=U(I) UV=U1(I) UAFU2(I) U(I)=OUT(I) U2(I)=SSO‘(U(I)‘UD)-SSZ‘UV-SS3‘UA U1(I)=UV+SSO‘UA§SS7‘U2(I) 896 CONTINUE C & C C TO DEFINE 111E DYNAMIC DEFLECTIONS IN mE LOCAL FRAMES C DO 100 I=2.17.3 J=I+1 IOO UL(I)3OUT(I)‘C3+OUT(J)'S3 C DO 200 I=3.18.3 J=I—1 2OO UL(I)=-OUT(J).S3+OUT(I)‘C3 C DO 300 1:21:3393 J=I+1 3OO UL(I)=OUT(I)‘C4+OUT(J).S4 C DO 400 1322:3493 131-1 400 UL(I)='OUT(J).S4*OUT(I)‘C4 C UL(36)=OUT(17)‘C4+OUT(18)‘S4 C RK11=(K-1)/1.0 C c: C THIS IS FOR THE PLO'ITING . IT C GIVES THE CRANK ANGLE. C UL(9) IS THE DEFLECTION IN THE MIDPOINT OF COUPLER C UL(28) IS THE DEFLECTION IN THE.MIDPOINT OF ROCKER C: WRTTE(7.153)RK11.UL(9).UL(28) (i 765 CONTINUE 153 FORMAT(F14.8.1X.F14.10) 113 CLOSE(7.STATUS='KEEP') STOP END 114 COOIOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.00.000.00.000. 00000 L O O L O O L O O L LLL OOOOO LLLLLLL THIS SUBROUTINE IS USED TO GET LINEAR LOCAL STIFFNESS AND MASS MATRICES 00000000000 El“??? cO...OO....00..OCO0.0BOOO...OOOOOOOOOOOOOOOOOOOOOOO SUBROUTINE LOL(E.A.EL.RI.MASS.KL.ML) REAL MASS.KL(6.6).ML(6.6).EL.E.RI C======= C INITIALIZE THE MY OF STIFFNESS MATRIX C======= DO 1 I=1.6 DO 1 I=1.6 1 KL(I.J)=0.0 c======= KL(1.1)=E‘A/EL KL(2.2)=12.‘E’RI/EL“3. KL(3.3)=4.‘E‘RI/EL KL(4.4)=KL(1.1) KL(5.5)=KL(2.2) KL(6.6)=KL(3.3) KL(1.4)='KL(1.1) KL(2.3)=5.'E‘RI/EL“2 KL(2.5)=-KL(2.2) KL(2.6)=KL(2.3) KL(3.6)=2.‘E‘RI/EL KL(5.6)=_KL(2.3) DO 2 I=2.6 Il=I-1 DO 2 J=1.Il.1 2 KL(I.J)=KL(J.I) C &====== C INITIALIZE mE MY OF MASS MATRIX C==a==s=== no 100 I=1.6 DO 100 I=1.6 IOO ML(I.J)=0.0 31:31:39.2: C TO CREATE THE LOCAL MASS MATRIX. c======= ML(1.1)=140. 103 115 ML(2.2)=156. ML(3.3)=4.‘EL“2 ML(4.4)=ML(1.1) ML(5.5)=ML(2.2) ML(6.6)=ML(3.3) ML(1.4)=70. ML(2.3)=22.‘EL ML(2.5)=54. ML(2.6)=-13.‘EL ML(3.5)=-ML(2.6) ML(3.6)=-3.‘EL“2 ML(5.6)=-ML(2.3) DO 103 I=2.6 Il=I-1 D0 103 J=1 .Il ML(I.J)=ML(J.I) DO 104 I=1.6 DO 104 I=1.6 104 ML(I.J)=ML(I.J)‘(MA88‘EL)/420. RETURN END 116 CO0.0000000000000000000000..O...OOOOOOOOOOOOOOOOOOOOOOOOt C O c GGGGG L 00000 L III N N ‘ C G L O O L I NN N ‘ C G GG L O O L I N N N ‘ C G G L O O L I N NN ‘ C GGGGG LLLLL OOOOO LLLLL III N N ‘ C O C THIS SUBROUTINE IS DESIGNED TO CONSTRUCTED ‘ C THE GLOBAL STIFFNESS AND MASS MATRICES ‘ C O COOOO.......‘COOOOOOC.....OO.............‘OOOOOUOOOCOOOOO C SUBROUTINE GLOLIN(X.XJ.X1.X2.X3.X4.XG.XM) REAL X(696)DXI(696)DXZ(6'6)Dx3(6'6)DX4(6D6)DXJ(636) REAL XG(40.40).XM(36.36) C C======: C INITIALIZE EVERY ENTRY OF THE MATRIX C======= DO 111 I=1.40 DO 111 J=1.40 111 XG(I.J)=0.0 C======= C CONSTRUCTED THE 40 BY 40 MATRIX WIHWT H.C. C====:== DO TO I=1.3 DO 10 I=1.6 IO XG(I.J)=X1(I.J) C DO 20 I=4.6 DO 20 J=4.6 K=I-3 M=J-3 20 XG(I.J)=X1(I.J)+X(K.M) C DO 30 I=4.6 DO 30 137.9 K=I-3 M=J-3 3O XG(I.J)=X(K.M) C XG(7.7)=X(4.4)+X(1.1) XG(7.8)=X(4.5)+X(1.2) XG(7R9)‘X(4s6)+X(103) C XG(8.7)=X(5.4)+X(2.1) XG(8.8)=X(5.5)+X(2.2) XG(8.9)=X(5.6)+X(2.3) C XG(9.7)=X(6.4)+X(3.1) XG(9.8)=X(6.5)+X(3.2) 16(9e9)=X(696)+X(303) 100 117 D0 40 I=7.9 DO 40 J=10.12 K=I-6 M=J-6 XG(I.J)=X(K.M) DO 60 I=10.12 DO 60 J=10.12 K=I-3 M=J-3 XG(I.J)=XG(K.M) DO 70 I-10.12 DO 70 J=13.15 K=I-9 MéJ-9 XG(I.J)=X(K.M) DO 80 I=13.15 DO 80 J=13.15 K=I-3 M=J-3 XG(I.J)=XG(K.M) DO 90 I=13.15 DO 90 J=16.18 K=I-12 M=J-12 XG(I.J)=X(K.M) XG(17.17)=X(5.5)+X2(2.2) XG(17.18)=X(5.6)+X2(2.3) XG(18.17)=X(6.5)+X2(3.2) DO 100 1816.18 DO 100 J819.21 K-I-15 M=J~15 XG(I.J)=X2(K.M) XG(19.19)=X2(4.4)+X3(4.4) XG(19.21)=X2(4.6) 110 1101 120 XG(20.19)=X2(5.4)+X3(5.4) XG(20.20)=X2(5.5)+X3(5.5) XG(20.21)=X2(5.6) XG(21.19)=X2(6.4) XG(21.20)=X2(6.5) XG(21.21)=X2(6.6) XG(19.37)=X3(4.1) XG(19.38)=X3(4.2) XG(19.40)=X3(4.6) XG(20.37)=X3(5.1) XG(20.38)=X3(5.2) XG(20.39)=X3(5.3) XG(20.40)=X3(5.6) D0 110 I=22.24 DO 110 I=22.24 K=I-21 M=J-21 XG(I.J)=X4(K.M) DO 1101 I=22.24 DO 1101 J=25.27 K=I-21 M=J-21 XG(I.J)=X4(K.M) XG(25.25)=X4(4.4)+XJ(1.1) XG(25.26)=X4(4.5)+XJ(1.2) XG(25.27)=X4(4.6)+XJ(1.3) XG(26.25)=X4(5.4)+XJ(2.1) XG(26.26)=X4(5.5)+XJ(2.2) XG(26.27)=X4(5.6)+XJ(2.3) XG(27.25)=X4(6.4)+XJ(3.1) XG(27.26)=X4(6.5)+XJ(3.2) XG(27.27)=X4(6.6)+XJ(3.3) DO 120 I=25.27 DO 120 J=28.30 K=I-24 M=J-24 XG(I.J)=XJ(K.M) XG(28.28)=XJ(1.1)+XJ(4.4) XG(28.29)=XJ(1.2)+XJ(4.5) 118 XG(28. 30)=XJ(1.3)+XJ(4.6) XG(29.28)=XJ(2.1)+XJ(5.4) 119 XG(29.29)=XJ(2.2)+XJ(5.5) XG(29.30)=XJ(2.3)+XJ(5.6) XG(30.28)=XJ(3.1)+XJ(6.4) XG(30.29)=XJ(3.2)+XJ(6.5) XG(30.30)=XJ(3.3)+XJ(6.6) DO 130 I=28.30 DO 130 J=31.33 K=I-27 M=J-27 130 XG(I.J)=XJ(K.M) DO 140 I=31.33 DO 140 J= 31.33 K=I-3 M=J-3 140 XG(I.J)=XG(K.M) DO 150 I=31.33 DO 150 J=34.36 K=I-30 M=J-30 150 XG(I.J)=XJ(K.M) DO 160 I=34.36 DO 160 J=34.36 K=I-3 M=J-3 160 XG(I.J)=XG(K.M) DO 170 I=34.36 DO 170 J=37.39 K=I-33 M=J-33 170 XG(I.J)=XJ(X.M) XG(37.37)=XJ(4.4)+X3(1.1) XG(37.38)=XJ(4.5)+X3(1.2) XG(37.39)=XJ(4.6)+X3(1.3) XG(37.40)=X3(1.6) XG(38.37)=XJ(5.4)+X3(2.1) XG(38.39)'XJ(5.6)+X3(2.3) XG(38.40)=X3(2.6) XG(39.37)=XJ(6.4)+X3(3.1) XG(39.38)‘XJ(6.5)+X3(3.2) XG(39.39)=XJ(6.6)+X3(3.3) XG(39.40)=X3(3.6) XG(40.40)=X3(6.6) 120 D0 180 I=1.40 DO 180 J=1.40 180 XG(J.I)=XG(I.J) 0:38:33 C PUT'B.C. IN THE MATRIX REDUCED TO 36 BY 36 MATRIX C:=a===s:a= DO 181 I=1.36 DO 181 151.36 181 XM(I.J')=O. C DO 190 I=1.19 DO 190 J=1.19 K=I+2 M=J+2 190 XM(I.J)=XG(K.M) C DO 200 I=20.36 DO 200 J=20.36 K=I+4 M=J+4 200 XM(I.J)=XG(K.M) C RETURN END 121 COO.........0.000.300.0000...OOOOOOOOOOOOOOOOOOOOOOOO 000000 0 O O O O O LLLLLLL OOOOOO F‘P‘P‘P‘ SZ:Z!Z!Z§§ 12 32i5522252 §§C>C’¢D§ C>C>C> EZEZIZSZii :Z 121E125232 THIS SUBROUTINE IS TO GENERATE NONLINEAR LOCAL STIFFNESS MATRIX 000.00.00.0000.00..0.0.0....OOOOOOOOOOOOOOOOOOOOOOOO KG(1-12)----NONLINEAR LOCAL STIFFNESS MATRIX C3fififififlflflfiflfififl C) SUBROUTINE LONON(E.A.UR.EL.UL.KG) REAL UR(6).KG(5.6).E C======= C INITIALIZE EVERY “TRY OF THE MATRIX C======= DO 830 I=1.6 DO 830 I=1.6 830 KG(I.J)=0.0 C KG(2.2)=6.0/5.0 KG(3.3)=2.0'EL“2/15.0 KG(5.5)=KG(2.2) KG(6.6)=KG(3.3) KG(2.3)=EL/10.0 KG(2.5)=-KG(2.2) KG(2.6)=KG(2.3) KG(3.2)=KG(2.3) KG(3.5)=-KG(2.3) KG(3.6)=’EL"2/30.0 KG(5.2)=_KG(2.2) XG(5.3)='KG(2.3) KG(5.6)=-KG(2.3) KG(6.3)=KG(3.6) RETURN END 122 COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.0.0.0...OOOOOOOOOOIOOOOOOOO (SCDC)CDCTCTCICICDCDCDCJCTCSCDCIC) 201 GGGGGG L OOOOOO NN N 000000 NN N G L O O N N N O O N N N G GGG L O O N N N O O N N N G G L O O N NN O O N NN GGGGGG LLLLLL 000000 N N 000000 N N BUILD UP GLCBAL MATRICESS (NONLINEAR PART OF THE PROGRAM) KM---GLOBAL MMTRICES FOR.EVERY ELEMENT KX--GL(BAL MATRICES FOR THE MECHANISM WITHOGJT B.C. KN---GLOBAL MATRICES FOR.THE MECHANISM WITH B. C. ......OOGOIOCOO OOOOOOOOO0.000000.000....0...000.00..OOOOOOOOOOOOOOOOO0.00......O. SUBROUTINE GLONON(KN.KM1.KM2.KM3.KM4.KM5.KM6. KM7.KM8.KM9.KM10.KM11.KM12.KX) REAL KN(36.36).KX(40.40) REAL KM1(6.6).KM2(6.6).KM3(6.6).KM4(6.6).KM5(6.6).KM6(6.6) REAL KM7(6.6).KM8(6.6).KM9(6.6).KM10(6.6).KM11(6.6).KM12(6.6) DO 111 I=1.40 DO 111 J=1.40 KX(I.J)=0. DO 10 I=1.3 DO 10 J=1.6 KX(I.J)=KM1(I.J) DO 20 I=4.6 DO 20 J=4.6 K=I-3 M=J-3 KX(I.J)=KM1(I.J)+KM2(K.M) DO 30 I=4.6 DO 30 J=7.9 K=I-3 M=J-3 KX(I.J)=KM2(K.M) DO 201 I=7.9 DO 201 J87.9 K=I-3 M=J-3 KK=I-6 MM=J-6 KX(I.J)=KM2(K.M)+KM3(KK.MM) 901 100 123 D0 40 I=7.9 DO 40 J=10.12 K=I-6 M=J-6 KX(I.J)=KM3(K.M) DO 60 I=10.12 D0 60 J=10.12 K=I-6 M=J-6 KK=I-9 MMBJ-9 KX(I.I)=KM3(K.M)+KM4(KK.MM) DO 70 I DO 70 J K=I-9 M=J-9 KX(I.J)=KM4(K.M) =10.12 =13.15 DO 80 I=13.15 DO 80 I=13.15 K=I-9 M=J-9 KK=I-12 MM=J-12 KX(I.J)=KM4(K.M)+KM5(KK.MM) DO 90 I=13.15 DO 90 J=16.18 K=I-12 M=J-12 KX(I.J)=KM5(K.M) DO 901 I=16.18 DO 901 J=16.18 K=I-12 M=J-12 KK=I-15 MMBJ-15 KX(I.J)=KM5(K.M)+KM6(KK.MM) DO 100 I=16.18 DO 100 J=19.21 K=I—15 M=J-15 KX(I.J)=KM6(K.M) KX(19.19)=KM6(4.4)+KM12(4.4) KX(19.20)=KM6(4.5)+KM12(4.5) KX(19.21)=KM6(4.6) 110 1101 120 1201 130 124 KX(20.19)=KM6(5.4)+KM12(5.4) KX(20.20)=KM6(5.5)+KM12(5.5) KX(20.21)=KM6(5.6) KX(21.19)=KM6(6.4) KX(21.20)=KM6(6.5) KX(21.21)=KM6(6.6) KX(19.37)=KM12(4.1) KX(19.38)=KM12(4.2) KX(19.39)=KM12(4.3) KX(19.40)=KM12(4.6) KX(20.37)=KM12(5.1) KX(20.38)=KM12(5.2) KX(20.39)=KM12(5.3) KX(20.40)=KM12(5.6) DO 110 I=22.24 DO 110 J=22.27 K=I-21 M=J-21 KX(I.J)=KM7(K.M) DO 1101 I=25.27 DO 1101 J=25.27 K=I—21 M=J-21 KK=I-24 MM=J-24 KX(I.J)=KM7(K.M)+KM8(KK.MM) DO 120 I=25.27 DO 120 J=28.30 K=I-24 M=J-24 KX(I.J)=KM8(K.M) DO 1201 I=28.30 DO 1201 J=28.30 K=I-24 MhJ-24 KK=I-27 MM=J-27 KX(I.J)=KM8(K.M)+KM9(KK.MM) DO 130 I=28.30 DO 130 J=31.33 K=I-27 M=J-27 KX(I.J)=KM9(K.M) DO 140 I=31.33 140 150 160 170 180 191 C 125 D0 140 J=31.33 K=I-27 M=J-27 KK=I-30 MM=J-30 KX(I.J)=KM9(K.M)+KM10(KK.MM) D0 150 I=31.33 DO 150 J=34.36 K=I-30 M=J-30 KX(I.J)=KM10(K.M) DO 160 I=34.36 DO 160 J=34.36 K=I-30 M=J‘30 KK=I-33 MMEJ-33 KX(I.J)=KM10(K.M)+KM11(KK.MM) DO 170 I=34.36 DO 170 J=37.39 K=I-33 M=J-33 KX(I.J)=KM11(K.M) KX(37.37)=KM11(4.4)+KM12(1.1) KX(37.38)=KM11(4.5)+KM12(1.2) KX(37.39)=KM11(4.6)+KM12(1.3) KX(37.40)=KM12(1.6) KX(38.37)=KM11(5.4)+KM12(2.1) KX(38.38)=KM11(5.5)+KM12(2.2) KX(38.39)=KM11(5.6)+KM12(2.3) KX(38.40)=KM12(2.6) KX(39.37)=KM11(6.4)+KM12(3.1) KX(39.38)=KM11(6.5)+KM12(3.2) KX(39.39)=KM11(6.6)+KM12(3.3) KX(39.40)=KM12(3.6) KX(40.40)=KM12(6.6) DO 180 I=1.40 DO 180 J=1.40 KX(J.I)=KX(I.J) D0 191 I=1.36 DO 191 I=1.36 KN(I.J)=O. COOOOOOOOOOOOOO C)C)C§C)f5 126 PUT BOUNDARY CONDITIONS IN OOOOOOOOOOOOOO DO 56 I=1.36 DO 56 I=1.36 IIM=I+2 JIM=J+2 IF(I.GE.20)IIM=I+4 IF (J.CE.20)JIM=J+4 KN(I.J)=KX(IIM.JIM) CONTINUE RETURN END 127 CO0.0000.00.00.00.0.00000000000000000000000000000000000000000000000 K K III N N K K I NN N K K I N N N K K I N N N K K III N N THIS SUBRWTINE IS DESIGNED TO CALCULATE THE ACCLERATION MD ANGULAR DISPLACEMENT L1.L2.13.IA GRGJND. CRANK. COUPLER AND ROCKER LENGTHES. (DCSCICDCICDCDCSCICDCSCDCTCDCICDCDCSCDCDCTCT SUBRWTINE KIN(RA. TH2 .TH3 .TH4 .TH21 . STEP. * XLEV1.XLEV2.AC.XLEN3.XLE€4.L1.L2.L3.IA) REAL L1 .L2.L3 .IA .AC(40) .RA(36.1) TH211=0.0 l A=2.‘L3"(L2‘COS(TH2)-L1) PI=3.1415927 B=L4"2-L1 ”2-L2 “2-L3 "2+2 ‘L1 ’L2 ’COS(TH2) C=2‘L2‘L3 ‘SIN(TH2) AZ=ABS(A“2+C'”2) D=SORT(AZ) AB(%ABS( (A‘BID“2 ) "2-(B"2-C“2) ID"2) DD=SORT(ABC) DA= 8IN(TH2) IF(DA.GT.0 . )mm TH3=ACOS(A‘BID“2+DD) ESE TH3 I=ACOS(A‘BID“2-DD) EIDIF m4=ACOS( (L2‘COS(TH2 )+L3‘COS(TH3 )-L1) ILA) TH31=-'m21"(L2‘SIN(TH2-TH4) )I(13 ‘8IN(TH3-TH4)) THETA2=TH29180JPI THETA3=m3‘180 . IPI THETA4=TH4‘180 . IPI TH41=-TH21" (L2‘SIN(TH2-TH3 ) ) /(L4‘SIN(TH3-TH4)) AA=TH21“2‘L2 *COS (TH2-m4) BB-TH3 1*‘2'L3 ‘COS(TH3-TH4) CC-L3 ‘SIN(TH3-'IH4) TH311=(TH31/TH21) ‘TH211-(AA+BB—TH41”2‘L4) ICC EE=TH21“2‘L2‘COS(TH2-TH3) FF=TH41 “2‘L4 ‘COS(TH3-TH4) PROGRAM IS BASED ON THE PAPER BY SMITH AND MAUNDER(1967) XLEN1.2.3.4 ARE THE LENG'mES OF JOINTS TH2 CRANK ANGLE IN RADS. WHILE TH21 CRANK ANG. VEOCITY IN RAD/S. ssssesesssssss C....QOQMIGQ OOOOOOOOOOOOOOOOOO0OOO.C.....OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 128 GG=L4‘SIN(TH3-TH4) TH411=(TH41ITH21)‘THle-(EE-FF+TH31*‘2‘L3)IGC C======= C AC CONTAINS NODAL ACCELERATIONS IN THE C GLOBAL FRAME. C======= AC(1)=-L2‘(TH21“2‘COS(TH2) + TH211“SIN(TH2)) AC(2)=-L2‘(TH21“2‘SIN(TH2)-TH211'COS(TH2)) C DO 1011 I=3.21.3 AC(I)=TH311 1011 CONTINUE C D0 1012 I=24.39.3 AC(I)=TH411 1012 CONTINUE C EL1=(L3-XLEN2-KLEN1)/4. DO 1013 I=7.16.3 II=(I-4) l3 J=I+1 AC(I)=AC(1)-(II‘EL1+KLEN1)‘(TH31“2‘COS(TH3)+TH311‘SIN(TH3)) AC(J)=AC(2)-(II‘EL1+XLEN1)‘(TH31“2‘SIN(TH3)-TH311'COS(TH3)) 1013 CONTINUE EL2=(L4-KLEN3-XLEN4)I4. D0 1014 I=28.37.3 II=(I-25)/3 J=I+1 AC(I)=~(XLEN4+II*EL2)‘(TH41“2‘COS(TH4)+TH411‘SIN(TH4)) AC(J)=~(XLEN4+II‘EL2)‘(TH41“2‘SIN(TH4)-TH411‘COS(TH4)) 1014 CONTINUE C AC(4)=AC(1)-XLEN1‘(TH31“2‘COS(TH3)+TH311‘SIN(TH3)) AC(5)=AC(2)-XLEN1‘(TH31“2‘SIN(TH3)‘TH311'COS(TH3)) C AC(19)=-L4‘(TE41"2‘COS(TH4)+TH411.SIN(TH4)) AC(20)=-L4‘(TH41“2‘SIN(TH4)-TH411‘COS(TH4)) C AC(22)=0. AC(23)=0. C C AC(25)='(XLEN4)‘(TH41“2‘COS(TH4)+TH411'SIN(TH4)) AC(26)=-(XLEN4)‘(TH41“2‘SIN(TH4)‘TH411‘C08‘TH4)) C AC(40)=TH411 C=s===== C TO INCREMENT mE CRANK ANGLE AND EENCE THE MECHANISM (#83333: 11 TH2=TH2+PII(180.‘STEP) C======= _ 129 RA IS DEFINED RELATIVE TO THE FINITE ELEMENT PRMRAM WHILE AC IS DEFINED REATIVE TO THE MEGANISM KINEMATICS PRWRAM.“ NW HAVE THE INTERFACE. O?OOOO DO 1015 I=1.19 J=I+2 RA(I,1)=AC(J) 1015 CONTINUE C DO 1016 I=20.36 J=I+4 RA(I.1)=AC(J) 1016 CONTINUE C RETURN DID 130 c...OOOOOOOOOOOO...OO0.00000000000000000000000000000 C a C S T R A I N ‘ C a C This snbrontin is designed to calculate ‘ c the strain of each element ‘ c e casesseeeeeeeeeseaseasseseeeeseeeeeeeaesesessseeeeee C SUBROUTINE STRAIN(HIS.BB.UL.EBSB) REAL EBS(1.1).BB(1.6).UL(6.1).EBSB(1.1) CALL MMLT(EES,BB.UL,1.6.1) RETURN END