FAILURE ANALYSIS OF BLAST LOADED DUCTILE BEAMSbyQmg LuoB.Sc., Dalian University ofTechnology, 1991A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIES(Department of Civil Engineering)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIASeptember 1994© QingLuo, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)_________________________Department of / V/ L &A;/7/A1The University of British ColumbiaVancouver, CanadaDate CA- / ‘ //IYDE-6 (2/88)IIAbstractA failure analysis procedure based on the elastic-plastic finite element analysis of thehigher order beam is presented to predict the mode II and ifi failure of impulsively loadedductile clamped beams.The variational equation of motion of the problem is first obtained by using the principleof virtual work and the Total Lagrangian Approach together with the kinematics of thehigher order beam theory and the nonlinear strain-displacement relation. The FiniteElement Method is employed to discretize the beam spatially initiating the numericalsolution procedure. The constitutive model considers elastic-plastic isotropic strainhardening by von Mises yield criterion and associated flow rule and strain rate sensitivityby Cowper-Symonds relationship. Equations of motion are integrated by centraldifference method in the time domain.Based on the finite element simulation, two failure criteria are proposed. In theinteraction failure criteria, contributions from tensile tearing and transverse shearing to thedamage of the beam are considered by including tensile strain ratio and shear stress ratio inthe criteria. Both linear and quadratic models are investigated. The beam is modelled asan elastic-plastic beam with a plastic hinge at each support to calculate the total strain atthe support. The shear stress comes either from the wall reaction obtained fromequilibrium or directly from the finite element analysis of the higher order beam theory. Inthe sectional plastic work density criterion, the beam is assumed not to fail until thesectional plastic work density exceeds the critical value. Post failure analysis is alsoincluded to take account of the residual energy of the beam at failure.111Numerical simulations have been carried out for the experiments of blast loaded beams.Comparison with the experimental results suggests the quadratic interaction criteria issuitable for this type of analysis.ivTable of ContentsAbstract iiTable of Contents ivList of Tables viList ofFigures viiAcknowledgments xi1. Introduction 11.1. Background 11.2. Purpose and Scope 62. Finite Element Formulation and Solution Procedure of the Problem 82.1. Introduction 82.2. Kinematics of the Higher Order Beam Theory 92.3. Variational Equation ofMotion 152.3.1. Variational Equation ofMotion 152.3.2. Governing Differential Equation and Boundary Conditions 192.4. Finite Element Formulation of the Equations ofMotion 212.5. Modelling the Material Behavior of the Beam 292.6. Solution Procedure Over the Time Domain 392.7. Summary 433. Failure Analysis of the Problem 453.1. Introduction 453.2. Interaction Failure Criterion 46V3.2.1. Calculation of the Influence of Tensile Tearing 463.2.2. Measure of the Influence of Transverse Shearing 543.2.3. Two Interaction Failure Criteria 583.2.4. Another Two Interaction Failure Criteria 593.3. Failure Criterion Using Plastic Work Density 603.4. Post Failure Analysis 643.5. Summary 654. Numerical Examples and Discussions 664.1. Introduction 664.2. Menkes and Opat’s Experiments 674.3. Modelling the Problem 694.4. Numerical Results and Discussion 724.4.1. Numerical Simulation of the Response Process 724.4.2. Discussion on Failure Conditions 804.4.3. Convergence and Energy Conservation 1134.5 Summary 1175. Conclusions 1215.1. Summary 1215.2. Areas Requiring Further Investigation 123References 125viList of Tables4.1 Convergence of mid-span permanent deformations with time step. 1164.2 Convergence of mid-span permanent deformations with mesh. 116vuList of Figures1.1 Three failure modes associated with the clamped metal beamsloaded impulsively. 32.1 Definitions ofpositive transverse displacement of the neutralsurface and rotation of a cross section of the beam at the neutralsurface. 132.2 Shear force edge effect for a cantilever beam loaded at x=0 witha tip load P. 222.3 An idealized bilinear model of the elastic-plastic material. 312.4 Idealization of elastic-plastic linear strain hardening materialbehavior with strain-rate sensitivity. 333.1 Typical deflection profile as determined by FENTAB. 483.2 Plastic hinge in a thuly clamped beam. 503.3 Shear distribution along the length of the half beam obtainedfromFENTABv.2. 554.1 Schematic drawing of experimental configuration. 684.2 Experiment results for series of 6061-T6 beams (0.25x1x8-in. beam). 704.3 Response profiles of one half ofa 0.375 in. by 1 in. by 8 in.beam subjected to an impulse of 52.9 ktaps. 734.4 Time history of central deformation for a 0.375 in. by i in. by 8in. beam subjected to an impulse of 52.9 ktaps. 754.5 Comparison of time histories of mid-span deformation obtainedfrom the FE analyses including shear (FENTAB v.2) and notincluding shear (FENTAB). 774.6 Comparison of time histories of mid-span deformation duringvmearly stage of response obtained from the FE analyses includingshear (FENTAB v.2) and not including shear (FENTAB). 784.7 Comparison ofbeam proffles during the early stage of responseobtained from FENTAB and FENTAB v.2. 794.8 Mid-span permanent deformations vs. impulses obtained from a0.375 in. by 1 in. by 8 in. beam using the quadratic interactionfailure criterion (QIC) and the linear interaction failure criterion (LIC). 814.9 Comparison of tensile strain ratios and shear stress ratios vs.impulses obtained from a 0.3 75 in. by 1 in. by 8 in. beamusing the quadratic interaction failure criterion (QIC). 834.10 Comparison of tensile strain ratios and shear stress ratios vs.impulses obtained from a 0.3 75 in. by 1 in. by 8 in. beam usingthe linear interaction criterion. 854.11 Comparison of pull-ins vs. impulses obtained from a 0.375 in. by1 in by 8 in. beam using the quadratic interaction criterion (QIC)and linear interaction criterion (LIC) with a ten element mesh. 864.12 Comparison of times of failure vs. impulses for a 0.375 in. by 1in by 8 in. beam using the quadratic interaction criterion (QIC)and linear interaction criterion (LIC) with a ten element mesh. 874.13 Mid-span permanent deformations vs. impulses obtained from a0.25 in. by 1 in. by 8 in. beam using the quadratic interactioncriterion and the linear interaction criterion with a ten element mesh. 894.14 Comparison of tensile strain ratios and shear stress ratios vs.impulses obtained from a 0.25 in. by 1 in. by 8 in. beam usingthe linear interaction criterion. 904.15 Comparison of tensile strain ratios and shear stress ratios vs.impulses obtained from a 0.25 in. by 1 in. by 8 in. beam usingixthe quadratic interaction criterion. 914.16 Time histories of total strain at the support for different finiteelement meshes. 924.17 Time histories ofwall reaction for different finite element meshes. 934.18 Comparison of time histories ofF and wall reaction for aneight element mesh. 964.19 Comparison of time histories ofF and wall reaction for antwenty element mesh. 974.20 Distributions of shear force along the halfbeam in the earlystage of response. 994.21 Mid-span permanent deformations vs. impulses obtained from a0.375 in. by 1 in. by 8 in. beam using the quadratic interactionfailure criterion with shear stress from the FE analysis (QICF)and the linear interaction criterion with shear stress from the FEanalysis (LICF). 1004.22 Comparison of tensile strain ratios and shear stress ratios vs.impulses for a 0.3 75 in. by 1 in. by 8 in. beam using the quadraticinteraction failure criterion with the shear stress from the FEanalysis (QICF). 1014.23 Distribution of shear force at 60 microseconds along half a 0.375in. by 1 in. by 8 in. beam subjected to an impulse of 17.8 ktapsfor a forty element mesh. 1034.24 Mid-span permanent deformations vs. impulses for a 0.25 in. by1 in. by 8 in. beam using the quadratic interaction failure criterionwith shear stress from the FE analysis (QICF) and the linearinteraction failure criterion with shear stress from the FE analysis (LICF). 1044.25 Comparison of tensile strain ratios and shear stress ratios vs.ximpulses for a 0.25 in. by 1 in. by 8 in. beam using the quadraticinteraction failure criterion with the shear stress from the FEanalysis (QICF). 1054.26 Mid-span permanent deformations vs. impulses for a 0.375 in.by 1 in. by 8 in. beam using failure conditions based onsectional plastic work density (SPWD). 1074.27 Distribution of sectional plastic work density along half a 0.3 75in. by 1 in. by 8 in. beam subjected to an impulse of 17.8 ktaps. 1094.28 Distribution of sectional plastic work density along half a 0.3 75in. by 1 in. by 8 in. beam subjected to an impulse of42.9 ktaps. 1104.29 Distribution of sectional plastic work density along half a 0.375in. by 1 in. by 8 in. beam subjected to an impulse of 64 ktaps. 1114.30 Time histories of different energies for haIfa 0.375 in. by 1 in. by8 in. beam subjected to impulses of(a) 17.8 ktaps (b) 42.9 ktaps(c)64ktaps. 1124.31 Kinetic energy for haifa broken beam of 0.375 in. by 1 in. by 8 in.subjected to different impulses. 1144.32 Kinetic energy for haifa broken beam of 0.25 in. by 1 in. by 8 in.subjected to different impulses. 1154.33 Time histories of energy difference ratio using different timesteps for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulseof 17.8 ktaps. 1184.34 Time histories of energy difference ratio using different timesteps for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulseof 52.9 ktaps. 119xiAcknowledgmentsThis work is possible due to the efforts of many people. Particularly, Professor Olsonhas given me advice, supervision, and guidance which have been invaluable. Also, I wouldlike to thank Dr. J. Jiang, Mr. Nagaraja Rudrapatna and Mr. Julien Fagnan for their helpfi.ildiscussions. Special thanks to Mr. Scott Tomlinson for his patient typing, proof reading,and moral support. Financial assistances were provided by the graduate researchassistantship from the Natural Sciences and Engineering Research Council grant and theUniversity ofBritish Columbia Graduate Fellowship.Finally, I am grateful to my parents for encouraging me to study in Canada. Withoutthem none of this would have been possible.1Chapter 1Introduction1.1 BackgroundIt is becoming increasingly necessary to examine the influence of impact loading onstructures in the design of various engineering systems. An understanding of the inelasticresponse of structures subjected to dynamic loads which cause permanent displacementsor structural damage is required for various branches of engineering. For example, toavoid the destructive action of earthquakes on buildings, such an understanding can beused to guide the development of rational design procedures. This knowledge is alsoessential to estimate the slamming and wave damage to ships and other marine vehicles.In addition it is applicable in studying the response of reentry vehicles and military vehiclessubjected to blast.Examining the influence of impact loading conditions on complete structures is quitecomplicated because the influence will depend on the type of structure on which theimpact load is applied. As a simplification, research is restricted to estimating theinfluence of impact loading on individual structural components, such as beams and plates.Refs. [1] through [5] give a detailed review in this area. The present work is concernedwith the inelastic response of beams subjected to impulsive loadings which causepermanent deformation or damage.2As noted in Ref. [1], in addition to producing an overall structural response of the beam,a dynamic load causes stress waves to travel through the thickness of the beam. When theimpact loads are severe enough this transverse propagation of stress waves can causefailure of the beam in the form of spalling. This type of failure occurs very early (usuallywithin microseconds of the initial impact) while the overall structural failure lasts longer bya few orders of magnitude. It is therefore customary to uncouple this early wavepropagation behaviour of the beam from the overall structural response. The presentwork focuses on the overall response of the beam by assuming the external dynamic loadis applied to the mid-surface of the beam instantaneously.Menkes and Opat [6] conducted a series of experiments on the overall inelasticstructural response of filly clamped metal beams with rectangular cross sections which aresubjected to transverse impulsive loadings. From the experiments, three failure modeswere identified as the load is monotonically increased, shown in Fig. 1.1. The failureswere characterized as: mode I, large permanent ductile deformations of the entire beam;mode II, tensile fracture of the extreme fibre at the supports; or mode ifi, transverse shearfailure at the supports. To date, extensive research has been done on blast loaded beams.However, most of this work deals just with large inelastic deformation, i.e. mode I failure,and only limited research has been done on the dynamic rupture of beams, i.e. mode II andifi failure.One of the few researchers who has studied the dynamic rupture of beams is N. Jones[7,8]. He proposed a simple rigid-plastic method of analysis, which considers theinfluence of finite displacements. This analysis, which accounts for the influence ofmembrane and bending forces, gave good agreement with the Menkes and Opat’sexperimental data from the beams which suffered mode I failure. To predict the onset of3(C)Fig 1.1 Three failure modes associated with the clamped metal beams loaded impulsively[8j. (a) mode I, large permanent ductile deformations. (b) mode II, tensile tearing atsupports. (c) mode ifi, transverse shear failure at supports.(a)(b)4mode II failure, the beam is modelled with plastic hinges at the supports and mid-span.After the maximum strain is calculated by using the model of plastic hinge, the thresholdcondition for a mode II failure can be obtained by equating the maximum strain to thestatic uniaxial rupture strain of the material. A square yield criterion relating themembrane force and the bending moment is used. To estimate the threshold for mode ififailure, the concept of transverse shear slide (which is analogous to the concept of aplastic bending hinge) is specified as an idealization of rapid changes of slopes across ashort length of the beam. The mode ifi failure is then assumed to occur when the amountof the transverse shear sliding at the supports reaches the beam thickness. Here, a squareyield condition relating the transverse shear force and bending moment is used. It is statedthat encouraging agreement with the experimental data for the onset of both mode II andifi has been obtained. However, Jones neglected elastic deformation and the effect oftransverse shear on deformation. The interaction effect between the tensile and shearaction (which both contribute to mode II and ifi failure) is also not included.In the rigid-plastic theory, as employed in Refs. [7,8], material elasticity is ignored whenthe external dynamic energy is significantly greater than the maximum amount of strainenergy which can be absorbed in a wholly elastic manner. However, according to Ref. [1],for material elasticity not to exercise a dominant influence, the external dynamic energymust be at least five times larger than the corresponding strain energy which can beabsorbed in a wholly elastic manner. Moreover, when the durations of external pressurepulses are comparable to the natural period of elastic vibration, the material elasticity canalso have a significant influence on the response of the beam.Fagnan [10] predicted the onset of mode II and III failures while accounting for theeffect of finite displacement based on the finite element procedure for the classical EulerBernoulli beam proposed by Folz [9]. To predict the mode II failure, Fagnan also used the5plastic hinge model for calculating the maximum strain at the support, and defines thatmode II failure occurs when the maximum strain reaches the rupture strain in the staticuniaxial tensile test. For the mode ifi failure, it is assumed that rupture of the beamoccurs when the shear stress exceeds the ultimate shear strength. Because analyses basedon the classical Euler-Bernouffi beam cannot estimate shear stress, it is calculated from thewall reaction force obtained from the equilibrium of the symmetric beam. Yet, this theorystill does not consider the transverse shear effect on the deformation of the beam and theinteraction between the tensile and shear action. Furthermore, the Euler-Bernoulli beamtheory is proved to be fundamentally inadequate for the analysis of impact conditions.Nonphysical results have been reported that waves of infinitesimal wavelength (whichmeans discontinuities) propagate with infinite velocity. It can be shown that satisfactoryresults are achieved only for lower modes and the shear-stress series for this theory arenon converging [6].According to Ref. [1], the transverse shear effect is believed to have a more importantinfluence on the response of beams when loaded dynamically than when loaded statically.Ref. [14] also points out that the effect of shear flexibility and rotary inertia play animportant role in theories of beam vibration and dynamic behaviour under impulsiveloading and it should not be ignored in the analysis. As to the interaction, no sharpdistinction is found between failure modes II and ifi. Failures involving both of thesemodes were observed when the beams were subjected to impulsive loadings in the rangebetween mode II and ifi thresholds. Therefore, interaction effects of both tensile andshear actions should be included in the failure conditions.A more universal energy criterion to predict the inelastic failure modes of beams loadedimpulsively and which includes effects of bending, membrane and shear action, issuggested by N. Jones and W.Q. Shen [11,12]. The rupture of the beam is assumed to6occur when the actual plastic work absorbed in a plastic hinge exceeds the critical value.The length of the plastic hinge is shown to vary with f3, the ratio of the plastic shear workto the total plastic work dissipated in the plastic hinge, according to an empirical relation.The mode ifi failure occurs when the ratio 13 reaches the critical value. As before, thismethod assumes a rigid-plastic material, which neglects the elastic deformation. Anempirical relationship between the plastic hinge length parameter and the energy ratio 13 isnot confirmed. Further experimental and theoretical studies are required to assess theaccuracy of the failure condition.The present research overcomes many of the above limitations. As explained later, thepresent work is based on elastic-plastic finite element analysis [13]. It considers the nonuniform transverse shear effect by using a higher order beam theory. And, it includes theinteraction effect of tensile and shear actions in the failure conditions.1.2 Purpose and ScopeAs stated in the previous section, the purpose of the present work is to develop a simplefailure criterion to predict the onsets of mode II and ifi failure for impulsively loadedbeams based on the elastic-plastic finite element analysis of the higher order beam [13].After the failure criterion is incorporated into the finite element analysis, responses(including failure mode) of beams loaded impulsively in the full range of intensities will besolved.Only the overall structural response of the beam is considered, which includes largeinelastic deformation and/or rupture. Within this research scope, initially straight7rectangular beams, with at least one axis of cross-sectional symmetry, undergo moderatelylarge planar deflection but small strains.The finite element formulation using the higher order beam theoiy is presented inChapter 2. To account for the non uniform shear distribution across the beam depth,higher order beam kinematics is used to describe the beam’s geometry. Since the stressstate is two dimensional, von Mises yield condition and associated flow rule as well asHooke’s law is used in the elastic-plastic constitutive relation. The Virtual Work Principleand Total Lagrangian Approach are utilized to obtain the variational equation of motion,and the finite element method is applied to determine the finite element formulation. Thespatially discretized equations of motion are then solved by integration over the timedomain with central difference. Chapter 3 discusses the modeling of the failure predictionof the problem. Two kinds of failure criteria are formulated. In the interaction failureconditions, interaction is exhibited in the way that the failure function is a mathematicalfunction of both tensile and shear contribution. In the failure criterion using the plasticwork density theory, interaction is considered in that the sectional plastic work densityincludes contributions from both tensile and shear stresses. To evaluate the accuracy andefficiency of the failure modeling, numerical simulations of the experiments conducted byMenkes and Opat are carried out. Results obtained from the simulations are compared tothe experimental results, and these are presented in Chapter 4. Finally, Chapter 5summarizes the research work, draws conclusions from the comparisons, and outlinesareas requiring further investigation.8Chapter 2Finite Element Formulation and Solution Procedure of the Problem2.1. IntroductionIn order to study the failure of the blast loaded beam, it is necessary to know thedisplacement field, stress, and strain state of the beam. This is accomplished by using theFinite Element Method. The solution algorithm, which is capable of numericallysimulating the transient dynamic response of the higher order beam undergoing geometricand material non-linear behaviour, is presented in this chapter. First, three different beamtheories are compared and the kinematics of higher order beam is studied in detail. Then,the Virtual Work Principle and Total Lagrangian Method are used to obtain the equationsof motion. After applying the Finite Element approximation of the displacement field, theFinite Element formulation of the equation of motion is obtained. The elastic-plasticmaterial model is also employed, taking account of plasticity by von Mises yield criterionand associated flow rule. Finally, the transient response of the beam is arrived at byintegrating the equations of motion using the central difference method.Although a similar formulation can be found in Ref. [9], the higher order beam insteadof the Euler-Bemouffi beam is employed in the present work to specify the beamkinematics. Also, the stress state here is of two dimensions instead of one, and von Mises’yield condition and associated flow rule are used to model the multidimensional stress-9strain relation in the present work. Moreover, the external damping is now taken intoaccount in the equation of motion.2.2 Kinematics of Higher Order Beam TheoryThe following solution procedure is only relevant to the bending effects of a slenderbeam undergoing large deformations. Large deformations are defined having themaximum lateral deflection experienced by the beam large relative to its height, althoughthe deflection is small relative to the longitudinal dimension of the beam. The torsion-freebending of the beam can be realized in the (x-z) plane by placing certain geometric andloading restrictions i.e. the cross-section of the beam has (x-z) plane as its longitudinalplane of symmetry and the resultant of transversely applied loads lies in this longitudinalplane of symmetry. For a rectangular beam of interest here, the x axis is set to coincidewith the centroidal axis of the beam and the z axis is chosen to be the transversesymmetric axis of the beam cross section. Since the longitudinal dimension of a slenderbeam is much larger than its lateral dimensions, it is assumed that the stress componentsa a, and can be neglected in comparison with the other stress components i.e.a=o=t=O (2.1)There are three theories to describe the beam kinematics. Euler-Bernoulli beam theoryassumes that cross sections normal to the undeformed neutral axis remain plane,undistorted, and normal to the deformed neutral axis. The shear deformation is thereforeneglected. Following this assumption, as well as those previously made, the displacementfield can be obtained10ôw(x)u1(x,y,z)= u(x)— zu2(x,y,z)=0 (2.2)u3 (x,y,z) = w(x)where u(x) and w(x) represent, respectively, the horizontal and transverse mid-surfacedisplacements of the beam. Shear deformation is accounted for in the Timoshenko beamtheory [19,20], which assumes planes which are normal to the beam axis in theundeformed state remain plane in the deformed state but no longer normal to the beam’sdeformed neutral axis. The displacement field obtained isu1 (x,y,z) = u(x) + z(x)2(x,y,z)=0 (2.3)u3(x,y,z)= w(x)where (x) is the rotation of the cross sectional plane. However, the relaxation of theclassical Euler-Bernoulli assumption in Timoshenko beam theory leads to thecontradiction that the resulting stress field no longer satisfies the shear-free boundarycondition on the surfaces of the beam. A correction factor, the Timoshenko shearcoefficient, is necessary to correct the contradictory shear stress distribution over the crosssection of the beam. To eliminate this deficiency of the Timoshenko theory, Levison [15],Bickford [16], and Heyliger and Reddy [17] introduced the higher order beam theory.This is a shear deformation theory for rectangular beams that retains the parabolicdistribution of the transverse shear strain, so there is no need to use the shear correctioncoefficients. The higher order beam theory correctly accounts for the stress free boundaryconditions on the upper and lower surfaces of the beam.As observed from the Timoshenko beam theory, the longitudinal displacement is only alinear expansion with respect to the thickness coordinate. In the current higher orderbeam theory, the displacement field is chosen to be of a more general form includinghigher order terms. Refs. [15] and [18] give a detailed description about this theory.11Because of the condition that the transverse shear stresses vanish on the upper and lowersurfaces of the beam and be non-zero elsewhere, the use of a displacement field in whichthe longitudinal displacement is expanded at least as cubic functions of the thicknesscoordinates is required. Since the beam is slender, both the Poisson’s ratio effects andstress components other than longitudinal normal stress and transverse shear stress areneglected. The transverse deflection of the displacement field can be set constant throughthe beam thickness. That is, the displacement field begins with,u1(x,y,z) = u(x) + zNJ(x)+z2(x)+z34(x)2(x,y,z)=O (2.4)u3(x,y,z)= w(x)where v is the rotation of a normal to the axis of the beam. The functions (x) and 4(x)are determined from the shear free condition on the top and bottom surfaces of the beami.e.t(x4)=O (2.5)This is equivalent to requiring the corresponding shear strain to vanish on these surfaces,giving1 &3 2 8W8=-j---+----=N!+2z3P+-- (2.6)Setting (x4) = o gives(x)=O and x)=—---(w+) (2.7)Therefore, the displacement field for the higher order beam theory becomes, if timedependence is taken into account,124z2 8wu1(x,y,z) = u(x,t) + z[W(x,t)— () (w+u2(x,y,z)=O (2.8)u3 (x,y,z) = w(x,t)This is with the positive transverse displacement of the neutral surface and the positiverotation of a cross-section of the beam at the neutral surface, shown in Fig. 2.1, and thedefinition of positive longitudinal displacement same as for the Euler-Bernoulli beamtheory. Indeed, 4(x), which is(2.9)can be called the warping function since it reflects the extent of the deviation of thedeformed beam cross section from an original plane surface. The displacement fieldimplies that, in the higher order theory, the cross section is not only allowed to rotaterelative to the neutral axis (as in the Timoshenko beam theory) but also to warp into anon-planar section, which is specified to satisfy the shear free boundary conditions at thetop and bottom of the beam.Since the beam undergoes finite deformation, Green’s strain tensor is employed todescribe the state of deformation of the beam. With the displacement field established, theassociated Green strain tensor (with reference to a Cartesian coordinate system) can beobtained by applying1 ãuj ôu. ôu ôu (i,j= 12,3) (2.10)to the higher order beam kinematics. Notice the restrictions imposed by assumptionsdiscussed previously. Non-zero Green strain tensor components follow as130__Fig. 2.1. Definitions of positive transverse displacement of the neutral surface androtation of a cross section of the beam at the neutral surface. From Ref. [15]‘(x,t) + -dxw(x,t)14=+!()2(2.11)£13 =(N’+)[1—4()2]Eq. (2.11) can also be put into this form:4z(z’26ii=e+zKi+--) 21(2.12)83 =[1_4()2]qÔU laW2where the axial stram e = - + (-i-) , the first curvature id = -, the second curvature6w= —(-k + —j-). and the effective shear stram ( = i41+ -.This is the non-linear strain-displacement relation for the higher order beam theory. Asobserved from the relation, the von Karman strain is included to describe the geometricnon-linearity caused by moderately large deflections of the beam. The effect of vonKarman strain is realized to be important when the transverse deflection is equal to orlarger than one thickness of the beam but still small relative to the length of the beam.Since the beam undergoes large deformation but small strains, the Green strain measure iseverywhere small, i.e.(2.13)152.3 Variational Equation of Motion2.3.1 Variational Equation of MotionWith the higher order beam kinematics specified previously, this section will establishthe variational equation of motion. Since the beam undergoes finite deformation and theGreen strain tensor is still small everywhere in the beam, a Total Lagrangian Approach isadopted. In this approach, the rectangular Cartesian coordinates are fixed in space as thereferences of the body before and after deformation, and the coordinates defining points ofthe original undeformed body are employed for locating the points of the subsequentlydeformed body. For this problem, the x axis of the fixed rectangular Cartesian coordinatesis set to coincide with the centroidal axis of the beam and the z axis is chosen to be thetransverse symmetric axis of the beam cross section in the undeformed configuration, asdiscussed in the previous section.The principle of virtual work is an avenue which leads to equations that govern thedynamic response of the structure. For the present dynamic problem with finitedeformation the virtual work principle can be written as [21]:IoV(s8u+°pu5u+°du6u)dv = SOS TöudS (2.14)This means that the work of external forces should be equal to the work of internal,inertial, and viscous forces given any small kinematically admissible motion. In the aboveexpression, °p,° lcd ,° V,° S represent, respectively, the mass density, the materialdamping parameter, the volume, and the surface area (on which the surface traction isspecified) of the beam before deformation. Sjj is the second Piola-Kirchoff stress tensor, uis the displacement vector, T is the prescribed surface traction vector, 6u is the small16virtual displacement vector which satisfies both compatibility and essential boundaryconditions, and the corresponding Green strain tensor. Each of the above terms isspecffied with reference to the undeformed configuration so as to be consistent with theTotal Lagrangian Approach concept. Here, the self weight of the beam has beenneglected because the inertia force in a dynamic problem is usually much greater.Following the previous kinematics discussion, for the higher order beam, the first term inthe virtual work expression (2.14) is.1 Su6EudV = .1(S1&+ 2S13&)dV (2.15)ov ovbecause s 1 and 613 are the only non-vanishing components of the Green strain tensor.According to Ref. [22], the second Piola-Kirchoff stress tensor is related to the Eulerianstress tensor Gjj (in which stress resultants are expressed) by the transformation equation°p ôu• ôuãu.S, =—[a — (i,j,k= 123) (2.16)where p is the mass density in the current deformed configuration, and ô is theKronnecker delta. For the present problem uj (i1,2,3) is set as for the higher order beamtheory discussed in the previous section (Eq. 2.8). Because the Green strain tensor iseverywhere small i.e. <<1, --- (1,k=1,2,3) is believed to be far less than 1 and°p p is considered to be valid. Therefore, S11 = a, S13 a3 = t, i.e.ISejdV= f(a&11 i-’&13)dV (2.17)ov ovIntroducing stress resultants,1771= SodA V= ft[1_4()2]dA?Th =fZ m;= S ()2d4j2(2.18)where 71 is the axial force, ml is the first order moment about the neutral axis, 771) is thethird order moment about the neutral axis, V is the effective shear force, t is the firstorder shear force, and t is the second order shear force. Obviously, V = V + t. Here, °Arepresents the original cross sectional area of the beam. Inserting the Green strain tensorcomponents and stress resultants into Eq. (2.17) gives= + fl-ö(-) +m18()—++ V(NJ+ .)jix(2.19)Eq. (2.19) can also be written as follows, using Eq. (2.12)Ssjoe, = f[noe + m1oic + + V&p]dx (2.20)where °L is the length of the undeformed beam.The second term on the left side of the virtual work expression Eq. (2.14), aftersubstituting the higher order beam kinematics (Eq. 2.8), is as follows [16,17]18J °püöudVDy= f(°puiöui+°pu3)dVOv168i’ 8w 18’ 8w=Jf p Auöu+j p INJ&V—j p I--NI—j p Iwô(--)-Fj p Iö(--)+pAw]dx°L(2.21)where 01 =-1—bh, the moment of inertia of the original beam cross section.12A similar expression can be written for the third damping term in Eq. (2.14).The work done by the external forces isITöudS = f[fou + fow]dx (2.22)Oswhere f and f are the distributed axial and transverse loads acting on the beam with unitsof force per unit length.Finally, summarizing the above discussion,J[no(.) + flö()+m1o() — + + Vö(j, +68 .. 16 8i 16 .. 8w 1 0’ 8w+° p°Aüöu p°Iöi.j,—j° p°I--ö—-j-° p°Iö(--) +-j° p°I--ö(--)+° p°A5w68 . 16 8w 16 . 8w1d°Aüöu j°1’d °IVöWjj°1Cd IöWj1cd °INJö()1 8*8wjd°Ijö(j)-1-°1d °A*6w—f8u — fow]dx =0 (2.23)This is the variational equation of motion for the higher order beam theory. Starting fromthis equation, a complete set of governing differential equations of motion, natural19boundary conditions, and essential boundary conditions can all be derived for the newbeam theory, which is briefly discussed next. However, more important is that from nowon the displacement field can be approximated by the finite element version and thealgebraic equation system for the unknown variables can be derived to solve the problemnumerically. This is significant as the analytical solution is almost impossible for thepresent problem, involving both geometric and material non-linearities.2.3.2 Governing Differential Equations and Boundary ConditionsWith the variational form of the equation of motion established, the governingdifferential equations will be persued in this section. This will be helpful when applyingboundary conditions to eliminate the rigid body motion in the Finite Element Method.After integrating the appropriate terms in Eq. (2.23) by parts and collecting thecoefficients of &i, ow, and 5j,, the governing differential equations are:= PAU+ICdAUam3 68 .. 16 ãi 68 . 16—-- — —i- + V = (—iöPIw+ jI--) + (—-jKdIw+ j’d’)o Ow O7l3 16 Ofr 1 O 16 0.j, 1(71 - + V) + f = pAM’ + p1- - p1 —i- + )CdAW +jd’E-(2.24)The appropriate boundary conditions are as follows:20natural essentialboundary boundaryconditions conditions71 u‘I’8w 8?712 16 1 8w 16 . 1 8Wfl————+V--—pI+—pI——-—1cI\+—1cjI— Wax 105 21 105 21 ax8wãxStress resultants are defined as in Eq. (2.18).As a consequence of generalizing the displacement field beyond the traditionalassumption (i.e. plane sections remain plane), the force resultants in the current higherorder theory are also more generalized. Third order bending moments and second ordershear forces occur which do not exist in the classical beam theory [161.According to the results of variational principle, since the force boundary condition isunknown at the clamped boundary, u w == N’= 0. Bickford [16] did threecombinations of boundary conditions to verify this. The same three combinations ofboundary conditions were checked in the present work: (1) u = w = =0, (2)u=w=V=0, and (3) u=w===0. Both results confirm thatu = w ===0 is appropnate for a clamped end. The reason for this is that since u,8w .w, and w on the centerline of the beam are used to descnbe the deformation of thewhole beam with the higher order beam kinematics (2.8), specifying uiO across the beamdepth at the clamped boundary makes it equivalent to forcing u = w == =0 tosatisfy ui=0 at every point on the clamped section.21For a simple supported boundary, from the variational principle, the essential boundarycondition should be uw=O. This is also confirmed by test calculations with combinationsof different boundary conditions.Also from the variational principle, the wall reaction force at a clamped boundary for astatic beam will be8w am3Rw=fljj+V (2.25)Actually this is also the transverse reaction force at the end of a static beam of anyboundary condition. If the von Karman effect is neglectedam3F=—---+V (2.26)is the appropriate transverse shear reaction force transmitted along the length of the beamfor this theory. This is verified by Bickford [16]. Furthermore, he found the boundarylayer character of the solution. For example, a cantilever beam is loaded with P at tip x=O.For this case, the boundary layer character of the elastic solution manifests itself in that atthe clamped boundary (x=L) V) approaches to 0 quickly, the higher order resultant m3increases rapidly to supply the total resistance, and the transverse shear force transmittedalong the length of the beam Fw remains equal to P, as shown in Fig. 2.2.2.4 Finite Element Formulation of the Equations of MotionHaving obtained the variational equation of motion, to solve the present probleminvolving both material and geometric non-linearities the Finite Element Method will beused to approximate the solution numerically. As usual the finite element approximationof the displacement field discretize the spatial domain of the beam into a number of22FPILFig. 2.2. Shear force edge effect for a cantilever beam loaded at x0 with a tip load P.From Ref.[161.23subdomains, or finite elements. Within each element the approximate displacement field isassumed in terms ofgeneralized nodal coordinates. According to the variational statementin Eq. (2.23), the transverse displacement should be at least twice differentiable and Ccontinuous, and the axial displacement u and rotation w should be at least oncedifferentiable and C° continuous [17]. Therefore, for a typical element, the displacementfield u is approximated by:(2.27)2 2 4u(q,t) = uj (t)Nj (rU, v(rI,t) = wj (t)N (q), w(1,t) = A, (t)M (TI) (2.28)Here, 11 is local axial coordinate, uj,j, Aj are nodal variables with A11,A3=w2 and for a typical element. Nj are the linear Lagrangian interpolationfunctions and Mj are the Hermitian cubic interpolation functions.Ni(rU=1—f N2(ri)=f= 131:11) +2(f) M2(TI)= 1i[1_2(f)+(f)] (2.29)= (J — 2(f) M4(fl) =where 4 is the length of the element. As can be seen, the displacement field u is a functionof both space and time, the shape functions are functions of space only, and nodalvariables are functions of time only. In this way the local separation of variables isaccomplished. In vectorial form, the appropriate set of generalized nodal coordinates foreach element is:1 ow iTdeUi,NJiWi() ,u21Vw1(J ). (2.30)X1 2J1e.Le___.1= 1ive,1e‘eue_Qttl ))—‘____________—£pUiv1eI1reNweiir 1eNeJ{o’o’o’11-’o}=n{io’o’‘io’o’-}=(t’€z)apt?=zia1a.+=dJ_ T[Z)Jj (EEc)Zi:si(j)uO!WnbuI14suJ9oui(j)U9W9jjoidAgjojpjjwowodsipUOTSJA1UUI1TUOtfl&irrnnsqnsoup1oopozqioujpoupajeToossJosuuupssaidxoquUUTtpiosuuiisuatp‘potsgqisoupjjuwJds!ppuTnssxuwwuowunjodqss![N]iAE(c)aP[N]=fl:sATwAuopUIT1OM11SITJS1JpunooZp11‘ivoo14t1(wz)ap[N]....apoo1N000TN0—Ih--fl000tN000TNL’J:swoaqiui.upuiojpjjuwo1Jds!putjj,qsowU!uwpjdAiqiimunbUATossa(.)uondusqnstj2500 0 0 00 0 000 0 0 00 0 00 0aM1 8M1 8M1 aM20 0ÔM1 M1 ôMiian ô’in a thiô0 06M2 ôM1 M2 aM20 0aM2 aM3 aM2 ôM4B-a1afl ai10 0 0 0 0 0 0 000 0 0 00 0 00 0ôM3 8M1 ôM3 ôM20 0aM3 ôM3 M3 8M4aa ao aa aa0 08M4 8M1 ôM4 8M20 08M4 ôM3 aM4 ÔM4aaq ao aaSince the beam is discretized into elements and the displacement field satisfies thenecessary continuity, the virtual work statement (2.14) for the discretized beams turnsinto:f (s,Jö6,+°pööu+icdüöu)dv= fTdudS (2.36)e=1 e=1where NE represents the total number of elements in the beam.For each element the work done by the internal forces:=J[n& + + moic2 + V&p]dii (2.37)where Ve aiid 6 are, respectively, the volume and the length of the element. With e, K1,1C2, and 4i defined previously in Eq. (2.34):I SjEj = 6deTI[flBi + ?mB2 + mB3 + VB + flBde]d (2.38)fSU&,ödeTrernt (2.39)Vewhere reul# = J[nB1 + m1B2 + mB3 + VB + flBde]dll (2.40)26is the internal force vector for a typical element which takes into account the internalresistance developed within the beam element in response to the external loading(regardless ofwhether or not the response is material and geometric linear). Its evaluationinvolves numerical integration techniques, which will be discussed later. Next isintroduced:retmL = + 7fl1B2 + m3B + VB4)dri re = IrLBdli (2.41)leThenreiW = retmL + (2.42)where reL is the linear part of the internal force vector of the element and rede is thenon-linear part. The presence of the non-linear part is because of the geometric nonlinearity.The work done by the inertia force in each elementlou TpiicIL?p°AüOu +-j p°NJONJ 1O5pOJfr5()+jp0IO() + pA17?OwIIx(2.43)After substitution of the finite element version of the displacement field, for a typicalelement in the finite element:IouTpudv= Ode Tmd (2.44)Vewhere me is the element consistent mass matrix. For the higher order beam theory it takesthe form of:27[mi’] [0] [0]me = [0] [m2] [m3] (2.45)[0] [m2] [m3]wherem1= IpAN1Ndri m=pIN1Ndn(2.46)23 C 16 dM• 32 33 fl_i dM1 dM 1me = J—p1d1Ndri =m efi me..=pdq dq+ pAM, M jdiiThe work done by the damping forces in each element can be similarly expressed ingeneral nodal coordinates in the above way, i.e.IuTKdudv 8dCede (2.47)VeThen the work done by the external force in a typical element can be written as:fouTfdq= 6dre (2.48)lewheref=jf,fI andret = 1[N1Tfdn (2.49)is the consistent external force vector for the beam element.Finally, in order to get a global set of generalized coordinates for the discretized beamfrom the element generalized coordinates, connectivity transformations in the form ofdeCed (2.50)are employed to illustrate the procedure. In fact, the Direct Stiffness Method is often usedin the practical application of the finite element method. Here d is the global vector of28generalized coordinates and C, is the connectivity matrix of bulean constants associatedwith each individual finite element under study.Therefore, the virtual work statement, as expressed by Eq. (2.36), turns out to be asfollows when expressed in the global generalized displacementsi3dT(CeTmeCed + CeTCeCe + CeTrt + CeTreCed) = ödT(CeTreet) (2.51)Since öd is arbitrary and kinematically admissible variations in the components of thegeneralized displacement vector, it can be cancelled out from the Eq. (2.51). So thegoverning differential equation of motion, as a result of the higher order beam theory andfinite element discretization, is.md+cd+r =r (2.52)NEwhere m = CeTmeCe is the global consistent mass matrix,e=1NEC = CeTCeCe is the global consistent damping matrix,e=1NErIt = C’r + CTrJ1/jCd is the global internal force vector, ande=1NErt = CeTret is the global consistent mass matrix,e=1The governing equation (2.52) can be interpreted as external loads are equilibrated by acombination of inertial, damping, and internal forces. The material does not have to beelastic since the stress-strain relation is not specified during the deduction.As observed, although the equations are spatially discretized, they are still a system ofsecond order ordinary differential equations in time, hence continuous functions of time.Therefore, they are called a finite element semidiscretization. An explicit direct integrationtechnique can be employed to discretize the system of equations in the time domain to29obtain a sequence of simultaneous algebraic equations [21]. This is discussed later inSection 2.5.2.5 Modeling the Material Behaviour of the BeamBefore starting the solution procedure for the governing differential equations of motionpresented in the previous section, rmt (the global internal force vector) must be evaluated.From the previous discussion, rmt can be obtained from remt (the internal force vector foreach element) by the Direct Stiffness Method or connectivity transformation. Eq. (2.40)can be used to derive rmt by using Gauss integration from stress resultants, i.e.fl, ?fl, ?fl and V, which again can be evaluated from Eq. (2.18) by using Gaussintegration if the stresses are known. So, to determine r1t, the appropriate constitutivelaw for the material must first be derived so that the stresses can be calculated from thestrains.IFor ductile materials subjected to high intensity transient loading conditions, plasticityand strain-rate sensitivity will inevitably be a consideration. It is therefore desirable toincorporate elastic-plastic strain hardening and strain rate sensitive material behaviour intothe material constitutive relation. In the current study a bilinear elastic-plastic strainhardening material model including the strain rate sensitivity is used for simplicity.Hooke’s law is used to describe the elastic part of the constitutive relation. Since thestress state is two dimensional, von Mises yield condition, the associated flow rule, andisotropic hardening are all employed to account for the plastic material behaviour of thebeam. Ref. [22] gives a detailed discussion on the theory of plasticity, which is veryhelpful.30In plasticity, the following basic assumptions are based on experiments:1. At any stage of deformation, the strain tensor can be expressed as the sum of theplastic strain tensor and elastic strain tensor , i.e.—(p) (e) 253where, 4, by its mean g, is related to the stress tensor by Hooke’s law,= (2.54)2. The plastic volume change is negligible. This is usually called the plasticincompressibility, i.e.=0 (2.55)To simpIi1r the material behaviour in the present problem, a bilinear elastic-plastic strainhardening model is used to describe the relation between stress and strain. Although thestress state in the present problem is no longer uniaxial, the bilinear model can still beapplied to the effective stress and strain. The theory of plasticity is later used to determinethe effective stress and strain. An idealized bilinear model of the material is shown in Fig.2.3.As observed from the bilinear model, the stress-strain curve is initially elastic (OA) withslope E yielding to a plastic strain hardening part (AB) with slope Et. That is, with astrain increment de, the stress increment is:da = Eds for aay (2.56)where ,as discussed before, is = + and de is the total strain increment, c1e isthe elastic component of the strain increment and is recoverable, de’ is the plasticcomponent and non-recoverable, and °y is the yield stress in uniaxial tests. The effect of31Fig. 2.3. An idealized bilinear model of the elastic-plastic material32strain rate on 0y will be addressed later. AlsoH= (2.57)Eis the plastic modulus. If Et is set to zero, H will be zero, resulting in an elastic-perfectlyplastic case. Upon unloading after yielding, the strain is reduced along an elastic path(BC), which is parallel to the elastic line OA; reloading retraces the unloading path, andfurther plastic deformation is produced when the maximum previously applied stress isexceeded.In addition to the bilinear elastic-plastic behaviour, the yield stress may now change inthe current problem. This is because the material is now subjected to high intensitydynamic loading, which causes high strain rates in the material. And it is shownexperimentally that, as the strain rate increases, the instantaneous yield stress of manymaterials increases from the quasi-static value. For some materials, this increase may bevery significant when the strain rate is high. To account for this effect, the CowperSymonds relation is used.=4i+(2.58)where °dy and Gy are, respectively, the dynamic and static yield stresses. D and p areexperimentally determined material strain rate parameters.Summarizing, the modeling of material behaviour is based on the ideal bilinear elastic-plastic linear strain hardening material models incorporated with strain rate sensitivity, asshown in the Fig. 2.4.Fig. 2.4. Idealization of elasticplastic linear strain hardening material behavior withstrain-rate sensitivity From Ref [23J33aa1a7eo=oe34To specify the constitutive relation, the yield condition must be known. Since the stateof stress is multiaxial in the present problem, a yield fi.inction has to be used to identifydifferent stages of material behaviour. Generally, the yield function depends on the stateof stress and strain and on the history of loading, i.e.41 (P),Kwhere ic, the work hardening parameter, is a function of the plastic strain tensor 41. VonMises yield condition is to be used in the current work. Introducing the equivalent oreffective stress as(2.60)with s representing the deviatoric stress tensor. In similar fashion, the equivalent oreffective plastic strain increment iscI = j2%d$d4 (2.61)The plastic work increment is dW = Since the material under study is isotropicwith strain hardening, the yield function can be written asf =f(o,H*()) (2.62)where H* ((P)) is the plastic stress-strain relation between & and Because vonMises yield condition incorporated with isotropic strain hardening is employed, f may bewritten asf =_H*()) (2.63)After yield, the way plastic deformation further proceeds is governed by the associatedflow rule, which states that the plastic strain increment must be proportional to thegradient of the yield function with respect to stress, i.e.35(2.64)Ii ôG,jwhere ? is the proportional coefficient. The hardening rule governs the way the yieldsurface (the f=O surface in the stress space) changes in size and shape as the plasticdeformation proceeds. For convenience, isotropic hardening (which assumes a uniformexpansion of the initial yield surface) is used in the current analysis, although it has littleexperimental support. The hardening rule for the present material model isH*(1)= b’ (2.65)which is specified by the bilinear model.As stated above, the yield thnction can be used to identifS’ different stages of materialbehaviour. With the yield condition well defined, the constitutive relation is ready to bedefined.For a state of stress corresponding tof(O, i.e. o< H, the material is elastic and there isno change in plastic deformation. The material stress-strain relationship is specified byHooke’s law which is, in vectorial form for the current two dimensional stress state, asfollowsdcr= [D]de (2.66)where [D] is the elastic constitutive matrix. This is applicable when the material is in theinitial elastic stage, it is unloaded, or is reloaded below the previous applied stress state.For a stress state with frO, i.e. i= H, the material yields and a change in plasticdeformation occurs. This is the case when the material is loaded beyond the yield state orwhen the yielded material, after unloading, is reloaded beyond the previously appliedstress state. Under these circumstances the stress-strain relationship is related to the36associated flow rule. Deduced from the theory of plasticity [23], the flow rule associatedwith von Mises yield condition in the form given in Eq. (2.63) isj8..C”) =(P)-i— (2.67)iiaGUIn vectorial form, the current two dimensional plastic strain state is:= (P) (2.68)where the effective plastic strain increment is[D]de=ãuj = (2.69)H+AA= 1—T {Df1—T (2.70)lôuJThus, whenf=O, de(P) and de(e) can be obtained from de with Eqs. (2.68), (2.69), and(2.70). By using Hooke’s law, the corresponding state of stress can be found. Theincremental stress-strain relationship may be briefly put asdo= [D]d6[D]tU [D]ID 1 = [D]— °1 (2.71)I H+Awhere [D,, j is the elastic-plastic constitutive matrix for the total elastic-plastic strainincrement.The state of stress wheref0, i.e. > H, has no meaning and is thus inadmissible. Tosummarize, Eq. (2.66) for yield functions less than 0 and Eq. (2.71) for yield functionsequal to 0 comprise the elastic-plastic constitutive relation.37Since the calculation of the equation of motion is carried out step by step (which isdiscussed later) sometimes after one time step of calculation f is greater than zero.Therefore, scaling is needed sincef>O is inadmissible. By scaling, stresses in the first statehave to be scaled to the yield state of stress, and the pure elastic portion of the strain hasto be subtracted, i.e.to=o2—a1 scal=—=todo-= —ds = 82—81 (2.72)= (oi) +scal*du (de) = (1— seal) *(d8)where (•)i and ()2 refer, respectively, to the given quantity associated to the first or thesecond state. After scaling, 01 is the stress state at yielding (which is a point on the yieldsurface) and de1 is the corresponding total strain increment (which is composed of bothelastic and plastic parts). Stresses can then be easily found using the plastic analysis forf=OinEq. (2.71).Having established the constitutive relationship, the state of stress at any point of thebeam may be determined given the strain. Attention now turns to solving for stressresultants i.e., 71, m1, m, V1. and V2 along the beam axis, which leads to the evaluation ofthe internal force vector. Because of the complicated distribution of these stress resultantsin elastic-plastic analysis, the integral expressions in which they are specified are bestintegrated numerically. Among the numerical integration schemes, Gauss quadrature isemployed in the current analysis because it is efficient, accurate, and suitable for therectangular beam cross section. When written with respect to the natural coordinate,=, the integration expressions for the stress resultants are as follows:38m1 = -Jc*i = (2.73)where b and h are the width and height of the beam respectively. Using Gaussianquadrature [241:‘==1T7f() (2.74)we have for the stress resultants:2 m 2= =i=1 1=1= 17t(j)= 1= V1 + t (2.75)which are enabled by the stress-strain relation presented before. Here, W1 is the weightassociated with the i-th depth-wise Gaussian evaluation point, and m=4. Four depth-wiseGauss points are sufficient to represent the non-linear stress distribution accurately as thematerial plastification extends through the beam cross section [9].After the stress resultants are obtained from the stress by Gauss quadrature, the elementinternal force vector can be easily evaluated from the integration equation (2.40) and alsovia Gauss integration. Written with respect to local coordinates, = — 1ret” = J(nB1+ + mB3+ + flBde)dC (2.76)using Gauss quadrature, gives:39reim = - wj[n(C)B1(C + ml(CJ)B2(Cf) + m(CJ)B3(Cf) + V(Cf)B4(Cf) + fl(Cj)B(Cj)de](2.77)where W is the weight factor associated with the jth spanwise Gaussian integration point.According to the definition of matrices, lBde requires the most number of Gauss points.Since the resultant axial force is going to approach a constant value under static load asthe finite element mesh is refined, three spanwise Gauss points are at least required withinan element in order to accurately represent flBde and hence relm. Therefore n=3 [9]. Withremt known, r” can easily be assembled by using the Direct Stiffness method.With all entries in the semidiscretized equation (2.52) known, the 2nd differentialequation system in the time domain may be solved, and this is discussed in the nextsection.2.6 Solution Procedure Over the Time DomainTo solve the dynamic differential equations over the time domain there are two types ofmethods: Modal Superposition and Direct Integration. For a blast loaded and non-linearproblem, the Direct Integration method is favoured. Among the Direct Integrationmethods, explicit time integration techniques are well suited to treat the material nonlinearities. This method can be implemented easily and can handle very large problemswith only modest computer storage requirements. The only disadvantage is that therequired time step to ensure computation stability is small [21]. However, in the presentblast loaded beam problem, the development of displacements and stresses are of interestand a small At is necessary for accuracy. Hence the explicit time integration is used in theanalysis. According to Ref. [9], the central difference is the most efficient method of time40integration when applied to problems involving impact or blast loading. Therefore, centraldifference is chosen to integrate the dynamic differential system of equations.In the Finite Element Method, the spatial domain is discretized into finite elements. As aresult, the finite element semidiscretization is:[mjfi}+[c]{á}+{rmt} = frt}By using central difference, the temporal domain of the problem is discretized into equaltime intervals of duration At. At each time interval a central difference approximation isused to replace the displacement time derivatives by differences of displacement {d):f} = — fd}_1) fJ} = ._--({d}÷1— 2fd} + {d}_1) (2.78)With the central difference method, the recurrence relationship [211 ism +—-c][d}÷ = {rt } — frmt} + ._!j..[m](2fd} — {d}_1)+(2.79)whereby the unknown displacements at t+At are given in terms of known quantities attime t. Since this is a two step method, a special starting procedure is required. Thestarting conditions are as follows:= fd}0 — &{a}0+.-fi}0= [m]_1({relt}o — {rh1t} — [c]fii}0) (2.80)As observed from the recurrence formulations, central difference clearly iscomputationally simple, since the unknown displacement vector occurring at any time step41is directly obtainable by performing simple arithmetic operations on matrices of knownvalue.If [m] and [c] are diagonal, then the equation system is uncoupled and the displacementvector can be obtained without solving simultaneous equations, thus making thecomputation even simpler. Actually, according to Ref. [21], explicit integration is usuallymore accurate with lumped mass matrices than consistent mass matrices. Therefore, withcentral difference used here for economy and accuracy, a lumped mass matrix is preferred.In the present study, the element lumped mass matrix is as follows:PA1y’io 0pAle3/r i_ /24LmeJ_ pAle//2340105PAl/pAle3//24(2.81)For the damping matrix, the extended Rayleigh damping matrix is used:[c]=a[m] (2.82)Therefore, both mass and damping are matrix diagonals, improving the economy andaccuracy of the computation.It is to be noted that central difference is conditionally stable and requires At such that:42AtAtcr (2.83)even when the structure is restricted to being linear. Here, Atcr is the critical time step toensure computational stability according to Courant’s criterionAtcr= 2 (2.84)0maxIf At is exceeded in time integration, the computations will be unstable resulting inerroneous time-histories as they grow unbounded. is the highest natural frequencyoccurring in the discretized linear or non-linear system. When the system is non-linear,because the geometric nonlinearity makes the system stiffer while the material nonlinearitydoes the opposite, an estimate of the critical time step is usually obtained by reducing thecorresponding linear result by 10-20% [9]. The present analysis can be further simplifiedbecause omax for the assembled finite element model is bounded by the maximum naturalfrequency of the smallest constituent unassembled and unsupported element, and thisfrequency can usually be computed easily [21]. For an Euler-Bernouffi beam, the criticaltime steps are as follows:for consistent mass matrix:Atcr = le/(q5cL) if le 1oq’Yr&cr = l/(10JJcLrG) otherwisefor diagonal mass matrixAtcrle/CL if le4.5rGItcr = l/(4...JcLrG) otherwise (2.85)where CL=is the speed at which longitudinal waves propagate along the beam axisVPrT. . .and rG=is the radius of gyration of the beam cross section. In the higher order beamtheory, shear deformation is taken into account and the beam is not as stiff as an EulerBernouffi beam, which means the highest natural frequency occurring in the higher order43beam will be lower. Therefore using the stability bound for Euler-Bernoulli beam elementis conservative, which is desirable.As observed from the above stability analysis, lumped mass matrix clearly provides forlarger stable time steps in addition to providing uncoupled equations and generally moreaccurate results than the consistent mass matrix in the central difference method. It is alsoto be noted that the critical time step in time integration is governed by the size of thesmallest element in the mesh, and the refining finite element mesh will require a morerestrictive critical time step. Therefore, extreme mesh refinement in any part of thedomain is undesirable.According to Ref. [21], the stability requirement is governing for the explicit timeintegration method. For systems of finite element equations, excellent accuracy can bearrived at by using a time step just under the stability limit. So in the present centraldifference method, a time step that satisfies stability criteria can usually guaranteeaccuracy satisfactorily.With the discussion about the integration scheme completed, the finite elementsemidiscretization may be turned into a system of uncoupled algebraic equations to solvefor the displacement vector. Moreover, by using the central difference, the solution in thetime domain may be used to obtain the dynamic transient response of the structure.2.7 SummaryThis chapter presented the formulation of the numerical scheme to simulate the blastloaded beam problem. Emphasis has been laid on the kinematics and the governing44equation of motion of the higher order beam theory, and elastic-plastic modeling of thematerial with strain rate sensitivity.With the whole formulation complete, the elastic-plastic transient response (includingthe non uniform shear effect) of the blast loaded beam may be numerically simulated. Thisprovides the basic data for the failure analysis later.45Chapter 3Failure Analysis of the Problem3.1 IntroductionUsing the formulation discussed in the previous chapter, the transient response of aslender ductile beam with geometric and material non-linearities under blast loadingconditions may be numerically simulated. With the displacement, stress and strain fields ofthe beam available, only the failure criteria is needed to predict failure when the beam issubjected to the fill range of load intensities. This is to be discussed next.In this chapter, several failure conditions are formulated. First, based on the previousresearch work by Fagnan [10], several versions of interaction failure criterion aredeveloped. Then, the failure condition using the theory of plastic work density is alsoformulated. Finally, the important post-failure behaviour of the beam is considered.All of the failure conditions are incorporated into the numerical simulation [13]discussed in the previous chapter. The results are presented in the next chapter fordiscussion.463.2 Interaction Failure CriterionMenkes and Opat [6] conducted a series of experiments on the dynamic plastic responseand failure of fi.illy damped metal beams which were subjected to explosive loadingconditions. At low impulsive loading conditions, large permanent ductile deformations ofthe beam were observed. At intermediate loading intensities tearing rupture wasencountered, which was caused mainly by excessive tensile strains developed at thesupports of an axially restrained beam. At still higher loading intensities, a more localizedfailure occurs owing primarily to the influence of large transverse shear forces. Betweenthe threshold values for the two failure modes, the beams failed with a mixed tensile-tearing and transverse-shear mode. This indicates that it is possible that both the tensile-tearing part and transverse-shear part contribute to the failure of blast loaded ductilebeams. This is the idea behind the interaction failure criteria. In this criterion, thecontribution of the tensile-tearing part is specified by a tensile strain ratio and the extent ofthe influence of the transverse shearing part is taken care of by a shear stress ratio. Bothparts are then included in the failure conditions to predict the failure of the beam.3.2.1 Calculation of Influence of Tensile TearingIn the interaction failure model, a tensile strain ratio is used to describe the contributionof the tensile-tearing behaviour of the beam to the failure of beam. Since the tearingaction of an axially restrained beam is most severe at the support, the maximum strain ofthe beam at the support is needed in the tensile strain ratio expression.47The elastic-plastic finite element analysis F.ENTAB version 2 [13,25], which isprogrammed according to the procedure presented in Chapter 2, is able to give thedisplacement, strain, and stress field of the present beam problem. A typical deformedprofile, obtained by FENTAB, of a clamped beam subjected to a uniform lateral load isshown in Fig. 3.1. However, as seen from the figure, the program is not able to accuratelymodel a plastic hinge at the support of the beam. A modified approximate shape of thebeam which accounts for the plastic hinge is therefore needed to predict the maximumtensile strain at the support.As mentioned in Ref. [26], Onat and Shield suggested that there exists a triangularplastic region at the support of the beam under impact conditions. This suggestion isencouraged by the appearance of deformed aluminum specimens in experiments [27].Thus, while using rigid plastic theory, Jones [7,8] modelled the clamped beam underimpulsive loading as two rigid parts joined by a plastic hinge at mid-span and two plastichinges at the supports. Based on this approach, Fagnan [10,27] derived a modelconsisting of plastic hinges only at the supports and the rest of the beam as modelled bythe Finite Element Analysis. This approximate shape of the beam is adopted in the presentwork when calculating the tensile strain at the support. In the elastic-plastic analysis,elastic-plastic deformation develops all along the beam. At the plastic hinge area at thesupport, it is assumed that the elastic response is negligible compared to the plasticdeformation.With the modified deflected shape of the beam established, the maximum total strain atthe supports can be developed from the beam shape model. In the present problem of afully clamped beam, the maximum total strain 6tot at the supports can be expressed asfollows6tot 8bend +8ial (3.1)48—— I—IL 9Fig.3. 1 Typical deflection profile as determined by FENTAB. From Ref.[1O]49where ej is the axial strain at the neutral axis at the support and 8nd is the bendingstrain at the support.Here, for a slender beam, the bending strain may be expressed asKh8bend = (3.2)where h is the depth of the beam and K is the curvature of the plastic hinge. For a plastichinge, the curvature K 15 defined as(3.3)‘phwhere 1’h is the length of the plastic hinge and 0 is the rotation of the plastic hinge.Because of the inability of the finite element analysis FENTAB version 2 to model aplastic hinge, approximation is made for the plastic hinge length 1ph A slip line fieldanalysis [26] can be used to obtain the shape of the plastic hinge across the depth of thebeam. It reveals that a plastic hinge in a beam has the shape shown in Fig. 3.2(a) when thebeam experiences pure bending behaviour and has infinitesimal transverse displacements.As the mid-span displacement increases, the membrane force increases and the axialstretching of the plastic hinge occurs until it takes the form in Fig. 3.2(b) with the onset ofa membrane response (when the mid-span transverse displacement equals the beamthickness h for rigid plastic analysis on a fully clamped beam). As observed from thefigure, it is evident that the length of the plastic zone on the upper surface of a fhllyclamped beam (iph) changes from beam depth h (when the beam is purely bending) to 2h(when the beam is at the onset of membrane state). Therefore it can be established that anapproximation for the hinge length “ph islph”th (3.4)50l=h(Q)lt2h____F -Fig. 3.2 Plastic hinge in a fblly clamped beam: (a) pure bending state; (b) onset ofmembrane state. From Ref.[12](b)51where a is the plastic hinge parameter whose value is chosen by the user. According tothe study done by Nonaka on the hinge length [26], for rigid plastic beams of rectangularcross section, 1 a 2 for maximum transverse displacement between 0 and h. After thatthe beam is regarded as being in the membrane state, the plastic deformation is assumed todevelop all over the beam, and the plastic hinge length is assumed to be:‘ph = whereL is the half length of the clamped beam. Jones [11,121 suggested that the plastic hingeparameter change with the impulse according to an empirical formula. In the presentelastic-plastic analysis, a is chosen to be 2 based on Fagnan’s research work [10].For the rotation of the plastic hinge, again an approximation has to be made because thefinite element procedure is not able to model a plastic hinge at the support. Theassumption is made that the rotation of the plastic hinge is equal to the maximum rotationoccurring in the beam, which is usually a short distance from the clamped support as seenfrom Fig. 3.1. To determine the maximum rotation of the beam, first the generalizednodal coordinates, as given by Eq. (2.30)I &w1 8w21rdeat all the finite element nodes can be examined, and the maximum nodal rotation valueoxcan be easily selected. Second, to be more accurate, interpolation of the nodal value canbe made along the length of each element using the shape functions to predict themaximum rotation. From Eq. (2.28) in Chapter 2.4 it is known that, within each element,the lateral displacement w has the relation to the element generalized nodal coordinates aswfO,0,Mi,M2,0, ,M3M4}de (3.5)whereM1,234are cubic Hermitian polynomials as specified in Section 2.4. Therotation of the beam can then be expressed by the first derivative of w with respect to x,52that is with respect to 11 (local element coordinate) within each element—10 0 8M1 8M2 0 0 8M3 8M4 ‘d 3 6Setting—i- = 0 and back substituting, the location and magnitude of the maximumrotation within each element can be determined. When written with respect to the localnormalized coordinate c 2’ —1, they are as follows.1(ãw‘2—_________________Cmax —6(w1_w2)+3le(-1+-)()max=— 1)W1 +(3C2 —2C— i) L+(i_c2)2+.4(3C2 +2C—(3.7)The interpolation procedure can be carried out for each element so that the maximumrotation for the whole beam, assumed to be 0, can be determined by comparison.Substituting the parameters discussed above into the bending strain formula (3.2) gives:0Ebend (3.8)2ciAttention now turns to the evaluation of the axial strain along the neutral axis at thesupport of the beam. This axial stretch can be obtained by directly integrating thedeformed profile as determined by finite element analysis of the beam element closest tothe support. Before deformation, the length of the element closest to the boundary is 10.During deformation, the arc length of the deformed element is then:5310 I1=0 ldx)After Taylor expression gives1 l0(f12 (3.10)llo +•Jfrom which the axial strain is obtained1 (3.11)8axial = 2lThis is the formula that is used to calculate axial strain at the boundary. Noticing the finiteelement approximation of the displacement field (3.5) within each element, for the firstelement:w=fO,0,M1,M2o,0 M34}d (3.12)where d1 is the displacement vector associated with the first element closest to theboundary. Then in the first element is the same asax100 8M1 8M2 00 ôM3 (3.13)Next110TEial_JdlT0,0 8M2 oM41 8M2 02l ( ‘o’aq’’’aq’aqJ(3.14)recognizingB — {o 0 8M1 8M2 0 aM3 aM4 1T 1 0 ôM1 8M2 0 0 8M3 ?jM4 1—Expressed with respect to the local nondimensionalized coordinate, C = —1, theICaxial strain at the support is_!‘6ax,al—J•diBdidc (3.15)54Choosing three spanwise Gauss points, which is required by the exact evaluation of B, theintegration by Gauss Quadrature is evaluated. In this way, the axial strain at the boundaryis determined.The total maximum strain at the support is now available by summing the bending andaxial components, as expressed in Eq. (3.1).Finally, the influence of the tensile tearing is measured by the ratio(3.16)where s is the nominal rupture strain of the material. This ratio can also be explained asthe contribution of the tensile tearing behaviour to the damage of the beam.3.2.2 Measure of the Influence of Transverse ShearMeasuring the influence of the transverse shear is accomplished by using a shear stressratio. According to the experiments, this influence is most significant at the support, andmetal beams are sheared off at the support during mode ifi failure. Therefore, tocompletely consider the influence of transverse shear, the shear stress at the support has tobe evaluated. In the present finite element analysis the higher order beam theory (whichconsiders non uniform shear distribution across the beam depth) is employed, allowinganalysis of the shear force and shear stress of the beam. As observed from the shear forcedistribution of the beam, however, the shear stress and resultants are more accurate in themiddle of the element than elsewhere. Shown in Fig. 3.3, the shear stress distributionfluctuates in the wave length of about one element size. This is probably because of theapproximate nature of the finite element method. To overcome this detriment of finite550 1 2 3 4 5Distance from the support (in.)2.0o i•1— Z01000800600.0C)0—2005Fig. 3.3. Shear force distribution along the length of the half beam obtained from FENTAB v.2.The beam is 0.375 in. by 1 in. by 8 in. Impulse is 17.8 ktaps.0 1 2 3 4Distance from the support (in.)56element analysis, two options are available. One is to obtain the shear stress by evaluatingthe reaction force at the support. The other is to use the shear stress at the middle of theboundary element to approximate the shear stress at the support. The first method shallbe discussed here. The alternative method will be dealt with later in Section 3.2.4.To get the shear stress from the shear reaction at the wall, the shear end reaction isrequired. This can be accomplished by applying the equilibrium in the vertical direction.Considering the beam as a free body, the vertical forces applied to the beam (including theexternal forces, inertia force, damping force, and end reactions) should be in equilibrium.Because of the symmetry of the beam, only half of the beam is considered, and theequilibrium can give the end reaction in this situation. The above discussion can be putinto the following expression.L L LR + jq(x)dx—Jm(x)dx—Jc(x)i4’dx= 0 (3.17)0 0 0where R is the wall reaction in the vertical direction, q(x) is the vertical component of theexternal distributed load, m(x) is the sectional mass density of the beam with a unit ofmass per unit length, c(x) is the sectional damping constant of the beam, and ‘ and *represent, respectively, the vertical acceleration and velocity of any point along the neutralaxis of the beam.As a result of applying the finite element method, the continuous beam is discretized intoan equivalent multi-degree of freedom system with the degrees of freedom {d}, associatedmass [ml, damping [c], and applied load {rt). Since lumped mass and damping areused, the equivalent multi-degree system reduces to a lumped mass system with mass anddamping associated to its own degree of freedom at each mode. Therefore, theequilibrium equation turns into:57R + q. — m131 — c13* =0 (3.18)nodes nodes nodeswhere qj is the effective external force on each node, mj3 and c13 are the mass anddamping associated with vertical degree of freedom at each node, respectively, and andw1 are the vertical acceleration and velocity at each node, respectively, which can beobtained from the transient finite element analysis. In this way, R can be found at eachtime step.With the shear reaction at the wall known, the shear stress is easily derived. For arectangular beam at an elastic stage, the maximum shear stress is at the centroid of thebeam section with a value of = 1.5Rw/ , where A is the area of the cross section.However, as the central fibers of the beam yields, the shear stress distribution across thedepth of the beam will become more uniform. In addition, it is assumed that the sheardominant failure will not occur until all the fibers in the beam cross section reach theultimate shear strength. Therefore, the shear stress in the present analysis is treated asbeing distributed uniformly across the beam depth at rupture of the beam, and can beobtained from the wall reaction as(3.19)Finally, to measure the influence of the transverse shear force, the normalized ratio(3.20)tultis employed, where tult is the ultimate shear strength of the beam when the beam is in thepure shear state. As with the tensile strain ratio, the transverse shear ratio reflects theamount of shear damage.583.2.3 Interaction Failure CriterionThe previous two sections discussed how to evaluate the influence of tensile tearing andtransverse shear on the failure of the beam. As stated ahead, a tensile strain ratio and ashear stress ratio are employed to measure the influences. These two ratios are alsoreflections of the damage extent caused by, respectively, tearing and shearing.With the measure of the two factors causing failure defined, the failure conditions maybe specified. The failure of the structure is related to tearing and shearing actions of thebeam. Therefore the failure condition must be a function of the two ratios, i.e.,(3.21)Ef tuft)Since both ratios represent the extent of damage to some degree, I ‘‘ can bejEj t)comprehended to be a general equivalent degree of the structure damage. As it isexpected failure will occur when the extent of damage is 100%, i.e. 1, it is then reasonableto assume that failure of the beam occurs when(3.22)Ef tft)For the failure thnctionf two functions are assumed:1) tmax = + tmaxI Ej t,jt ) Sf tuft(3.23)‘‘2,2)Sf taft) 8.1) taft59which are the linear interaction failure criterion (LIC) and the quadratic interaction failurecriterion (QIC).3.2.4 Another Two Interaction Failure CriterionAs mentioned is Section 3.2.2, the shear stress is obtainable from the finite elementanalysis because the shear effect is included in the analysis procedure. However, the shearstress is not accurate at the support because of the approximate nature of the finiteelement method. Since sufficiently accurate shear stresses at the mid point of the elementare obtainable, it is possible to approximate the shear stress at the support using the shearstress at the middle of the boundary element. In this way, as the finite element mesh nearthe support is refined, the value of shear stress at the wall is approached.The shear stress distribution across the beam depth is more uniform when the centerfibers in the section yield. It is also assumed the shear dominant failure does not occuruntil all the fibers in the section reach the ultimate shear strength. In addition, there mightbe possible fluctuations in shear distribution across the beam depth in the results of finiteelement analysis. Thus the average stress is usedJtcL4tavg (3.24)where V=ftdA, the first order shear resultant force.ASubstituting tavg for in the two failure conditions (3.23) in the previous section,another two interaction failure criteria are derived601) f1L, tayg = +1 Cf tuft) 6f tjjf(3.25)I 2,.. 22) +‘‘j Cf taft) I 8f) taft3.3 Failure Criterion Using Plastic Work Density TheoryThe work density caused by the stresses at each point of a structure consists of twoparts: density of strain energy, stored in the structure itself and density of plastic work,dissipated in the form of heat, etc. It is assumed in Ref. [11] that when the density ofplastic work reaches a critical value, a crack will initiate at that point. That is when(3.26)where & is the density of plastic work including the contributions from all the stresscomponents, and cr is the critical value of plastic work density. It can be taken as4?)= f Od(C,6)dC’ (3.27)where ad(C,é) is the dynamic engineering stress, which is obtained from a dynamicuniaxial tensile test for a given é, and where s is the plastic strain at rupture. Sinceelastic strain is far less than plastic strain at failure, plastic strain at failure is taken to beequal to the rupture strain, i.e.4?)=(3.28)where e is the engineering rupture strain. Actually ad and both could depend on themagnitude of the strain rate, in general. However, Cf is assumed to be independent of ê inthe present analysis. The Cowper-Symonds relation is used to take strain rate sensitivity61into account in a (E,). Since a bilinear material is assumed for the elastic-plastic strainhardening material, 0cr can be expressed as follows:cr dyef +--Es1. (3.29)where Et is the tangent modulus of the material and aê can be obtained from the CowperSymonds relation discussed in Chapter 2(flYa,=a 1+1—‘ DAs to the density of plastic work, S can be similarly expressed asE(p)= fajcic”” (3.30)where ad and e(P) are true dynamic stress and plastic strain, respectively, in an uniaxialcase. In an actual structure, they are assumed to be equal to the equivalent stress andstrain. In the present finite element analysis the stress and strain state at Gauss points areavailable at each time step. Thus, 0d and zSE can be obtained at each time step. S, thedensity ofplastic work at a point, can therefore be evaluated by trapezoidal integration.8(p)= (3.31)where aj1 and a1 are the stresses at i-1-th and i-th time step, respectively, and zSz is theplastic strain increment at i-th time step.When the density of plastic work is greater than the critical value, a crack is initiated atthat point. However, after a crack is initiated but before severance occurs, the structural62member might continue to support loads. The beam is therefore assumed not to fail untilthe whole section fails, that is, until the plastic work density in the whole section reachesthe critical value.cl=.clcr (3.32)where 2 is the sectional plastic work density, whose unit is work per length. cr is thecritical value.Because the plastic work density at each Gauss point is known, the sectional plasticwork density can be calculated. The sectional plastic work density is the plastic workdensity in the whole section, i.e.c=JJ&L4 (3.33)AFor a rectangular beam with width b and height h,f1=b8dz (3.34)Expressed in terms of local nondimensionalized coordinates,bh 1 (3.35)where = . Then the sectional plastic work density can be evaluated by ussquadrature. Four Gauss points are used in the Gauss integration. Therefore,bh ‘ (3.36)2 i1where 9. is the plastic work density at i-th Gauss point. W() is the weight factorassociated with i-th Gauss point.63Similarly, the critical value for the sectional plastic work density can be written ascr =JJcr’4 (3.37)AWhen expressed in terms of local nondimensional coordinates,bh 1cr=fcr’1 (3.38)2 1Because at each point of the section the strain rate might be different, the dynamic yieldstress at each point might be different and so is the critical plastic work density. Thereforeto evaluate ncr’ Gauss integration is employed. Applying Gauss quadrature,bh ‘2cr (3.39)2where is critical value for plastic work density at i-th Gauss point.Finally, with both sectional plastic work density 2 and the critical value cr established,the failure criterion using the sectional plastic work density can be stated as the beam willfail when, at one section,(3.40)This comparison between 2 and cr can be carried out at each spanwise Gauss pointalong the beam to decide whether the beam fails or not.643.4 Post Failure AnalysisWith the finite element procedure established in Chapter 2, and the failure conditionsdiscussed in the present chapter, the transient response and failure of the blast loadedbeam can be simulated. However, it should be noted that the transverse displacement ofthe beam at the instant of rupture is usually different from the permanent transversedisplacement obtained from the experiments. This is because the beam has a residualenergy at severance when all of the initial energy input has not been absorbed in plasticdeformations completely before failure, which is usually the case. The residual energy atrupture is often significant. With so much energy at failure, the beam is bound to deformuntil a stable state of the beam is reached. As the beam continues to deform afterseverance, part of the residual energy is absorbed in fhrther plastic flow until it reaches thesteady state where no more plastic deformation occurs. At the steady state, the beammoves with associated rigid body kinetic energy and vibrates with relatively small amountsof energy exchange between elastic energy and kinetic energy, simultaneously. The mid-span transverse deformation at the steady state, which is the average of the maximum andminimum value during the vibration, can be obtained from the post failure analysis.The idea of the post failure analysis involves treating the broken beam as a free-freebeam with initial motion. After severance of the beam at failure, the beam is blown offand becomes a free flying beam. To numerically simulate this, the boundary which wasclamped before is now set free. Further temporal progression of finite element analysis fora free-free beam can give information on the motion.After incorporating the post failure analysis into the program, it is now possible tosimulate the whole process of the motion of a beam subjected to blast loading conditions,from deformation to rupture to motion of a free beam. The comparison with experimental65data obtained by Menkes and Opat[6] and the discussion of the numerical model arepresented in Chapter 4.3.5 SummaryFailure conditions are developed in this chapter. The interaction failure criteria and thefailure criterion using plastic work density theory both consider all contributions fromtensile and shear action of the beam.In addition, to fully simulate the response of a blast loaded beam, and to check with theexperimental results properly, post failure analysis is discussed.This concludes the theoretical modelling of the problem. With the theoretical basepresented in Chapters 2 and 3, the whole response process of the beam may now besimulated numerically.66Chapter 4Numerical Example & Discussions4.1 IntroductionWith the finite element procedure discussed in Chapter 2, and failure conditions and postfailure analysis established in Chapter 3, it is possible to obtain the transient response ofthe blast loaded beam, to predict the failure, and to obtain the mid-span transversedeformation in the steady state. However, the efficiency and accuracy of the prediction isstill unknown.Menkes and Opat [6] conducted a series of experiments on 606 1-T6 clamped aluminumbeams subjected to a surface explosive charge. These experiments are used as the basisfor comparison. To check the efficiency and accuracy of the numerical model, finiteelement numerical simulations of these experiments are carried out. Calculation results arethen compared to the experimental data, giving encouraging correlations. Furthermore,some insights into the response of blast loaded beams are obtained through the finiteelement analysis.674.2 ExperimentsMenices and Opat [6] conducted a series of experiments on clamped aluminum beamssubjected to blast loads in which the uniformly distributed impulsive load is provided bysheet explosives. The experimental configuration is shown in Fig. 4.1. As shown in theconfiguration, the sheet explosive is cemented to a neoprene buffer which is bonded to thetop surface of the beams. Although single-ended detonation was used in the experiments,no significant evidence of deformation asymmetry is present in the experiment results, anduniformly distributed impulse loads can be assumed. The intensity of the impulse wasvaried by using different combinations of sheet explosive of different thickness. The beamspecimens are all of 6061-T6 aluminum. Three thicknesses (0.187, 0.25, and 0.375 in.)and two lengths (4.0 and 8.0 in.) were used. All beams were 1.0 in. wide.For the beams examined, as the load is monotonically increased, three distinctly differentdamage modes can be identified for the beams with rectangular cross sections, as shown inFig 1.1:I Large inelastic deformation of the entire beam,II Tearing (tensile failure) of the beam material at the extreme fibers at thesupports, andifi Transverse shear failure of the beam material at the supports.Within the load range of mode I, the beams responded in a ductile manner and acquiredexcessive permanent transverse displacements without any material failure. The damageseverity can be measured by the residual central deflection. At mode II the beams faileddue to dominantly tearing of the beam material at the supports. The threshold is taken asthe impulse intensity which first causes tearing of the beam material at the support.During mode II failure as the impulse increases the permanent deformation of the severedcentral section decreases (which can only be measured from photographs) until mode ifi is68— EXPLOSIVE— NEOPRENEBEAM______iI2LFig. 4.1. Schematic drawing of experimental configuration. From [6,10]69reached. This is shown in Fig. 4.2. In mode ifi, failure of the beam occurs because theinfluence of the transverse shear forces dominates the response and causes failure of thematerial at the support. The threshold is taken as the impulse intensity which first causesno significant deformation in the severed central section. Plastic deformation of the beambecomes localized near the supports and no appreciable plastic deformation is present overmost of the beam at this stage. The permanent deformation of the severed central sectionremains almost constant. When the beam is subjected to loads which lie between themode II and mode ifi thresholds, failures which involve both the tearing and shearingmodes are observed.4.3 Modelling of the ProblemThe modified finite element analysis FENTAB is used to numerically simulate the blastloaded beam experiments conducted by Menkes and Opat . For simplicity, the actualexplosive loading condition can be modelled as a rectangular pressure pulse, with aconstant peak value Po for a short duration t, distributed over the entire span of the beam.The impulse of the pressure pulse should be the same as the experimental value. Thevalidity of this simplification is because when Po’Pc is very large and the impulse is high,where Pc is the static plastic collapse load per unit length, the maximum transversedisplacement is insensitive to this dynamic pressure ratio. This suggests that the details ofthe pressure pulse are not important [1,81. According to Fagnan [10], the duration of therectangular pressure pulse is chosen to be 60 sec.Since the load and the clamped beam are symmetric, only a half span of the beam isstudied. To represent the half span in the calculation, several finite element meshes havebeen used. Ten elements of equal lengths are used most often. When a finer grid is70KIAPS10.917.825.628.735.139.642.946.050.752.96 .6Fig.4.2. Experiment results for series of 6061-T6 beams (0.25x1x8-in. beam) From [6]71needed, four different meshes were used: (i) sixteen elements of same size, (ii) twentyelements of equal length, (iii) twelve elements with the first two elements of the ten equalelement mesh subdivided into four equal length elements, and (iv) fourteen elements withthe first two elements of the twelve element mesh fi.irther subdivided into four equal lengthelements.The material of the beam is modelled as elastic perfectly plastic, i.e. Et=0. The reasonfor the modelling is that (according to Jones [11]) during the dynamic response of thebeam the strain rate will decrease with time, resulting in a stress-strain curve whichapproaches perfect plasticity. Other properties of the 6061 T6 aluminum as quoted byFagnan [101 are p=2.56x10 lbs2/in4,0o41,600 psi, and E=10.5x106psi. Since 6061T6 aluminum is essentially strain rate insensitive at the strain rates encountered in theexperiments, i.e.=o,. In the current study, damping is set to zero because the mostsevere condition of deformation is desirable for the analysis.To predict the failure of the beam, the ultimate data of the material is needed. Themaximum allowable strain is chosen to be 0.17, which is the percentage elongation atfracture of a static uniaxial tensile specimen made of 6061 T6 aluminum. Menkes andOpat [6] observed from the experiments that the dynamic shear strength appears to bemuch higher than the static one. After considering the von Mises yield criterion andultimate tensile strength of the material [10], the ultimate shear strength is chosen to be30,000 psi.With the experimental conditions modelled, results can be obtained from the numericalsimulation, which are discussed in the next section.724.4 Numerical Results and Discussion4.4.1 Numerical Simulation of the Response ProcessAfter modelling the problem of a blast loaded beam as described in Section 4.3, usingthe finite element procedure which is presented in Chapter 2, incorporated with failureconditions discussed in Chapter 3, it is now possible to numerically simulate the wholeresponse process of the beam subjected to blast loading conditions. An example of this isdiscussed below.Fig. 4.3 shows a series of proffles of one half of a blast loaded beam. The beam is 0.375in. high, 1 in. wide, and 8 in. long. The impulse is 52.9 kLaps (1 ktap = 0.014 lbf-s/in2 or1 ktap = 100 N-s/rn2). As observed from Fig. 4.3, the part of the beam close to theboundary undergoes severe deformation while the middle part of the beam experienceslittle deflection in the early response. As time progresses, the deformation graduallyextends to the middle of the beam. When certain failure conditions are satisfied, the beambreaks away from the supports, becoming a free beam. While flying away from thesupports, the beam experiences further plastic deformation, resulting in more energy beingabsorbed by the plastic deformation. Finally, the beam reaches a steady state of motion,and the beam moves with a rigid body motion superimposed by a vibration mode.The time history of the central deformation Aw of the same beam is shown in Fig. 4.4.The central deformation Aw is the mid-span deflection before rupture of the beam, and isthe difference of the displacement of the mid-span node and the boundary node afterrupture. This time history shows the whole response process from a different perspective.At rupture, the central deformation Aw is only 0.2 12 in. After the beam reaches thehighest peak, more energy is absorbed by plastic work, and then the beam reaches the73CC0E000I— 02C=0.150.0: 0.10C00.053:532.51.50.500.050Distace from origin (in.) 3 4— 0.04CCI00.020C00.0100250Dista,ce from origin (in.)00 1 2 .. 3 4Distance from origin (in.)Fig. 4.3. Response profiles of one half of a 0.375 in. by 1 in. by 8 in. beam subjected to animpulse of 52.9 ktaps.740 1DistLice from origin Qn.)0.70.6C— 0.5C0E00.400.30I:0CC0E00.0000C0ICC0E00a0.0000>0,CaI-4543.5 —32.5 -2u=0 40659 in.0550 microseconds, close to 1St highest peak of central deformation1Distnce from origin fin.) 587.5 —7—6.5 -8- u=0.24904 in. 1050 microseconds, close to 1st lowest peak of central deformation5.5 I0 1 .2.. 3 4Distance from ongin fin.)Fig. 4.3. Response profiles of one halt of a 0.375 in. by 1 in. by 8 in. beam subjected to animpulse of 52.9 ktaps. (continued)52.52CCo 1.5.(310.5075Fig. 4.4. Time history of central deformation for a 0.375 in. by 1 in. by 8 in. beam subjected toan impulse of 52.9 ktaps. The failure is predicted by the quadratic interaction failure criterion(QIC).0 1 2 3 4 5 6 7 8 9 10 11Time (millisec.)76steady state. The permanent deflection, which is taken to be the average of the maximumand minimum central deformation Aw in the steady state, is 1.474 in. Obviously, afterrupture, the beam experiences further severe deformation because of residual energy at therupture of the beam, which cannot be neglected. In this case, because of rupture of thebeam, the permanent deflection is less than the deflection when the beam is considered tohave indefinite tensile strain.Further study on the steady state of motion reveals that the period of the free beambroken away from the supports, as obtained from numerical results, is 820 .tsec. For auniform Bernouffi-Euler beam with both ends free, the circular frequency of first non rigidbody mode is [28]:____CD2 = (4.73 1)2.Vi- (4.1)which turns out to be 820 .tsec. Since the mid-span deflection time history of the beamwithout rupture obtained from the higher order analysis is almost the same as from theBernoulli-Euler beam analysis, as shown in Fig. 4.5, there should be no significantdifference between the natural periods of the two beams. The natural frequency of ahigher order beam with both ends free should be around 820 p.sec, which agrees with theresult from the numerical simulation.Compared to analysis without considering the shear effect, the time history of mid-spandeformation before rupture is almost the same, with one of the time histories shown in Fig.4.6. However, there are some differences between the profiles of the beams in the earlystage of response. Shown in Fig. 4.7 are some proffles during the early stage of responseof one half of a 0.3 75 in. x 1 in. x 8 in. beam subjected to an impulse of 42.9 ktaps. Sincethe shear effect is taken into account, the higher order beam is more flexible and therefore(as observed from the profiles) the part of the beam close to the boundary deforms more0.80.70.6CE0C00.20.10Results from FENTAB —-- Results from FENTAB v.2770.7Fig. 4.5. Comparison of time histories of midspan deformation obtained from the FE analysesincluding shear (FENTAB v.2) and not including shear (FENTAB). The beam is 0.375 in. by 1 in.by 8 in. and subjected an impulse of 17.8 ktaps.0 0.1 0.2 0.3 0.4 0.5 0.6Time (millisec.)0.50.40.20Co0.10__ Results from FENTAB _— Results from FENTAB v.2780.09Fig. 4.6. Comparison of time histories of midspan deformation during early stage of responseobtained from the FE analyses including shear (FENTAB v.2) and not including shear (FENTAB).The beam is 0.375 in. by 1 in. by 8 in. and subjected an impulse of 42.9 ktaps.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08Time (millisec.)1.50.50.040.03CC:0.02C)0=000.0100.250.20.15C0000.100.050079..... Results from FENTAB _..,.. Results from FENTAB v.2_.. Results from FENTAB ... Results from FENTAB v.2Dista3ioe from origin (in.)... Results from FENTAB .. Results from FENTAB v.2Fig. 4.7. Comparison of beam profiles during the early stage of response obtained fromFENTAB and FENTAB v.2.80severely in the early response. After the different responses in the early time, thedifference of the profiles close to the boundary between the two analyses decreases as timeprogresses. Finally, if the beam can sustain the impulsive load, the difference will be verysmall.4.4.2 Discussion of Failure ConditionsFrom the above discussion, the response processes of a beam subjected to differentloading conditions are now available. However, different failure conditions give differentpredictions of time of failure, or even different predictions ofwhether the beam can sustainthe impulse or not. That is, different failure criteria predicts different processes ofresponse. Therefore the accuracy of different failure conditions in predicting the failureneeds to be addressed. This is discussed in this section.1. Interaction Failure CriteriaAfter analyzing the 0.375 in. x 1 in. x 8 in. beam subjected to different impulsive loadsby using the interaction failure criteria, the mid-span permanent deformations (which arethe differences between the displacements of the mid-span and the supports) are plottedagainst the impulses in Fig. 4.8. Experimental results of mid-span permanent deformationis available from Ref. [6] for the beams which did not fail. Results for the broken beamsare digitized from the figure in Shen and Jones [11], which they obtained from measuringthe photographs in Menkes and Opat’s experimental report. For comparison, these resultsare also included in the figure.As observed from Fig. 4.8, the trend of the curve of mid-span permanent deformation,predicted by linear interaction failure criteria is correct. First, the mid-span permanentdeformation increases as the impulse increases. Then, after a certain impulse value, the2.5100Fig. 4.8. Midspan permanent deformations vs. impulses obtained from a 0.375 in. by I in. by 8in. beam using the quadratic interaction failure criterion (QIC) and the linear interaction failurecriterion (LIC).81C0(U10 20 30 40 50 60 70Impulse (ktaps)QIC (10 elements)QIC (20 elements)± Mode I Exp. ResultsLIC (10 elements)QIC (16 elements)x Shen & Jones’ (Mode II Results)82mid-span permanent deformation starts to decrease until the beam reaches mode ififailure, where the deformation is almost constant. However, this criterion predicts modeII tensile rupture of the beam far too early at 21.8 ktaps, as opposed to 46 ktaps reportedby Menkes and Opat. It also underestimates the mid-span permanent deformation of thebroken beams. Therefore, the linear interaction failure criterion cannot predict theresponse of the blast loaded beam as accurately.Using the quadratic interaction failure condition, results of failure analyses employingthree different finite element meshes (10, 16 and 20 equal elements) are presented in thefigure. The curves also predict the trend of mid-span deformation correctly. It fits theexperimental data very well and is on the conservative side. A detailed look at the failureof the beams reveals that the first rupture of the beam occurs at 289 IJ.sec. According toMenkes and Opat, confirmed by analyses of the beams which sustained the impulsive load,the maximum values of central deflection are reached at about 300 sec for the 8 in. beam.Hence, the failure analysis predicts the time of first rupture very reasonably. However, thequadratic interaction failure criteria still gives a fairly early mode II rupture. Experimentaldata ofMenkes and Opat shows that the mode II tensile rupture occurs at 46 ktaps for the0.3 75 in. x 1 in. x 8 in. beam, while the quadratic interaction failure criteria predicts themode II threshold to be 28.7 ktaps.To predict the mode ifi failure, the tensile strain ratios and the shear stressratios are plotted versus impulses. The results of three quadratic interaction failureanalyses employing different grids are presented in Fig. 4.9. As seen from the figure, thetensile strain ratios are greater than the shear stress ratios and they are almost constantduring the mode II failure. Then, at a certain impulse, significant increases or decreasesappear in the curves of tensile strain ratio or the shear stress ratios. The ratios of shear8310.9 tensile strain ratio0.8U,00.7 -0.6shear stress ratio0.50.4I I I I20 30 40 50 60 70Impulse (ktaps)_._ tensile ratio (10) _-•-_ tensile ratio (16)—A-- tensile ratio (20)__- shear ratio (10) __ shear ratio (16) __ shear ratio (20)Fig. 4.9. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained from a0.375 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion (QIC).84stress exceed the ratios of tensile strain and start to be the controlling factor in the failureof the blast loaded beam. The impulse immediately after the big increase or decrease isspecified to be the threshold of mode III. For the 0.375 in x 1 in x 8 in. beam, it is 50.7ktaps, while the threshold for mode ifi failure is reported to be 64 ktaps. Thisdiscrepancy is probably because of different definitions of mode ifi threshold. The modeifi threshold in the present analysis means a starting point where the shear stress ratiostarts to dominate, while the mode ifi threshold in Menkes and Opat’s experiment is takenas the impulse intensity which first causes pure shear failure with no significantdeformation of the severed central section. The strain and stress ratio chart for analysisusing absolute failure criterion is shown in Fig. 4.10, and 39.6 ktaps is given as the modeifi threshold, which is far too early.Using a 10 equal element model of the beam, the pull-ins and times of failure of the blastloaded beam obtained from the two analyses incorporating different interaction failurecriteria are plotted versus the impulses in Figs. 4.11 and 4.12, respectively. As shown inFig. 4.11, the pull-in versus impulse curves exhibit similarities to the curves of mid-spanpermanent deformation versus impulse. This is intuitive because the beam will be pulled inmore as the beam deforms more severely laterally, provided the beam is inextensible or haslimited extensibility. In Fig. 4.12 the first predicted time of failure for both curves occursaround 300 sec, which is reasonable as explained previously. During mode ifi failure,the predicted times of failure are very early in both curves (less than 60 p5cc, which is theduration of the pressure pulse). This also agrees very well with the experimentalobservation that mode ifi failures occur early. However, it is clear that the times of failurepredicted by the linear interaction criterion are far earlier than the ones predicted by thequadratic interaction failure criterion.8520 30 40 50Impulse (ktaps)Fig. 4.10. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained from a0.375 in. by 1 in. by 8 in. beam using the linear interaction criterion.0.90.8 -0.7 -0.6 -0.4 -0.3 -0.20.1 I I I I_— tensile strain ratio (10 elem)60-_.._ shear stress ratio (10 elem)70860.70.60.5. 0.4C0.30.20.10._. QIC (10 elemets) __ LIC (10 elements)Fig. 4.11. Comparison of pull—ins vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beamusing the quadratic interaction criterion (QIC) and the linear interaction criterion (UC) with a tenelement mesh.0 10 20 30 40 50 60 70Impact (ktaps)87350300a 250a)Cl)2C)2004l50 -100500-0 10 20 30 40 50 70Impulse (ktaps)_. QIC (10 elements) .... LIC (10 elements)Fig. 4.12. Comparison of times of failure vs. impulses for a 0.375 in. by 1 in. by 8 in. beam usingthe quadratic interaction criterion (QIC) and the linear interaction criterion (LIC) with a tenelement mesh.LIC6088The same failure analysis procedure incorporated with two different interaction failurecriteria has also been carried out for a 0.25 in. x 1 in. x 8 in. beam subjected to differentimpulse loads. The mid-span permanent deformation predicted by the two analyses usinga 10 equal element mesh are compared with the experimental data in Figs. 4.13 and 4.14.As observed from the figure, both curves predict the trend of mid-span permanentdeformation very well. However, the linear interaction failure criterion gives the mode IIthreshold as early as 17.8 ktaps (from Fig. 4.13) and mode ifi threshold as early as 28.7ktaps (from Fig. 4.14). According to the experimental data, 28.7 ktaps is reported to bethe mode II threshold and 46 ktaps the mode ifi threshold. Comparatively, the quadraticinteraction failure criterion gives a better prediction with mode II threshold as 21.8 ktaps(though this is still early) and the mode ifi threshold as 35.1 ktaps, as determined from thetensile strain ratio and shear stress ratio shown in Fig. 4.15.To check the mesh effect on the failure condition, different grids are used in calculatingthe response for the case of a 0.375 in. x 1 in. x 8 in. beam subjected to an impulse of 17.8ktaps, which (from the experiments) is supposed to cause mode I failure. The timehistories of the maximum total strain at the support and wall reactions are shown in Figs.4.16 and 4.17.In Fig. 4.16, the time histories of the maximum strain at the support obtained by fourdifferent meshes, i.e. 8, 10, 16, and 20 equal element grids, are presented. The four timehistories are observed to be very similar, and the differences between the curves decreasesas the mesh becomes denser. That is, the maximum strain at the support is independent ofthe mesh and converges as the mesh is refined. Therefore, the tensile influence portion ofthe failure condition, which is described by tensile strain ratio will converge withrefining of the mesh.2.5CC0EIa)4-.Ca)CCua)0C.O.502QIC (10 elem)A Experimental Results-_-_ LIC (10 elem)x Shen & Jones’8970Fig. 4.13. Midspan permanent deformations vs. impulses obtained from a 0.25 in. by 1 in. by 8in. beam using the quadratic interaction criterion and the linear interaction criterion with a tenelement mesh.0 .10 20 30 40 50 60Impulse (ktaps)10.90.80.70.60.4 -0.3 -0.2 -0.1010—— tensile strain ratio (10 elem) .÷ shear stress ratio (10 elem)9070Fig. 4.14. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained froma 0.25 in. by 1 in. by 8 in. beam using the linear interaction criterion.20 30 40 50 60Impulse (ktaps)91001.110.90.80.70.60.50.40.30.220 30tensile strain ratio (10 elem) . shear stress ratio (10 elem)70Fig. 4.15. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained froma 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction criterion.40 50 60Impulse (ktaps)920.090.080.07t00.06Cl)0. 0.050.040.0100 0.1—.-— 8 elements0.2 0.3 0.4 0.5 0.6 0.7Time (millisec.)_-_- 10 elements —i--- 16 elements0.8 10.9—-- 20 elements1.1Fig. 4.16. Time histories of total strain at the support for different finite element meshes.936543C 0o o2cOCu zQ)0h—2—30 0.1——8 elements0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9Time (millisec.)—— 10 elements —-- 16 elements —— 20 elements1.1Fig. 4.17. Time histories of wall reaction for different finite element meshes.94Comparisons of wall reactions calculated from different finite element meshes arepresented in Fig. 4.17. The same four meshes as used above are employed. During mostof the time history, the four curves are almost the same. At the very early stage of theresponse (0 - 75 sec) there is a small difference, but as the mesh is refined the resultsconverge. The results from the 16 element grid are very close to the results calculatedfrom the 20 element grid. When the stable response starts (500 j.i.sec) there is a small shiftbetween the responses calculated from the different meshes. This is probably because theprincipal frequencies of the different finite element lumped mass models of the beam areslightly different. Yet, as shown in the figure, the time history of the wall reactionobtained from 16 element mesh is almost exactly the same as the one from the 20 elementmesh. That is, with refining of the mesh, the results converge. From the point of view ofthe failure criterion, it is the wall reaction force in the early state of response which playsan important role. As previously discussed, the results of wall reaction force areindependent of the element size and converge when the mesh is refined. Therefore, theshear stress ratio is independent of mesh size, as is the shear portion of the/ tulifailure criterion.Summarizing the above discussion, both tensile and shear parts of the failure criterionare independent of the mesh and converge as the mesh becomes denser, and therefore theprediction obtained from the interaction criteria are independent of meshes and converge.The above detailed study at maximum strain at the support and the wall reaction forcehelp to explain Fig. 4.9. This figure, during the second failure mode, shows the tensilestrain ratio and shear stress ratio as almost the same regardless of the mesh. This isbecause the maximum strain at the support, which is dominant, is independent of elementsize and the wall reaction force is almost the same for different meshes between 75 - 50095sec. When the beam starts to fail in mode Ill, the ratio curves start to be different. Thisis because although the wall reaction force (the controffing factor in mode ifi failure) isalso mesh independent, there are differences in results calculated from different meshes inthe early time history (0 - 75 i.tsec), which is when the mode ifi failures occur. From Fig.4.17, it is known that the wail reaction force converges as the mesh density increases.Therefore, as Fig. 4.9 shows, the curves of the tensile strain ratio and shear stress ratioconverge as the mesh is refined during the mode ifi failure part.2. Another Two Interaction Failure CriteriaAs discussed in Chapter 2, from variational principle, the relation between the resultantforces and wall reaction force is obtained.6W ã?73 16 .. 1 6W 16 . 1 6W= 71—— + V-—pI’iJ + —pI—-—lCdIW (4.2)Thus, in a case when the wail reaction force is not available from the equilibrium (such asan asymmetrically supported beam) it can still be calculated from the resultant forces.This relation is checked for in the case of the 0.375 in. x 1 in. x 8 in. beam subjected to animpulse of 17.8 ktaps. In this case, as stress related results are more accurate in middle ofthe element, the resultant forces at the middle of the boundary element are used tocalculate F, the transverse shear resultant force transmitted along the length of the beam,to approximate the wall reaction R.(4.3)6x axBecause the midpoint of the boundary element is quite close to the boundary, theacceleration and velocity related terms are neglected in the calculation. As the elementmesh is denser, the F should be able to approach the wail reaction. Comparisons ofF(the transverse shear force) and the wall reaction force are presented in Fig. 4.18 for an 8element mesh and Fig. 4.19 for a 20 element mesh. As observed from the graphs, for both966543.0Cl)w o2—2—31.10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9Time (millisec.)Fw —--- Wall reactionFig. 4.18. Comparison of time histories of Fw and wall reaction for an eight element mesh.9754-3-20 D0)h -—23 I I I0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1Time (millisec.)Fw__Wall reactionFig. 4.19. Comparison of time histories of Fw and wall reaclion for a twenty element mesh.98meshes the Fs at the mid point of the first element are very close to the wall reactionforces. For a denser mesh (e.g. 20 element mesh) the approximation is slightly better. Aspreviously discussed, the wall reaction will converge as the mesh is refined. Hence theresultant force, including the shear stress, will converge as the refining element mesh.Shown in Fig. 4.20 is the distribution of shear forces in the middle of each element alongthe beam at seven different times. The beam and the loading condition are the same asmentioned above. A 20 element mesh is used to model the beam. It can be clearlyobserved how the shear forces develop during the loading. When the loading stops at 60p.sec, the shear force at the boundary reaches its maximum.Using the interaction failure criterion with the shear stress calculated from the finiteelement analysis, analyzing the 0.3 75 in. x 1 in. x 8 in. beam subjected to the whole rangeof impulsive loads results in the mid-span permanent deformation versus impulse figure asshown in Fig. 4.21. As before, the linear interaction failure criteria predicts the mode IIthreshold earlier than the quadratic interaction one, and underestimates the mid-spanpermanent deformation during mode II failure. For the quadratic interaction failurecriterion, results from three meshes are compared in the graph. All three curves fit thedata well and predict the increasing-decreasing-constant trend of mid-span permanentdeformation. The mid-span permanent deformation improves when the mesh is refined to12 elements from 10 elements. This is because stresses obtained from the program areused in the analysis and a finer mesh is needed to give accurate enough results, as opposedto an analysis based on displacement related quantities. The mode II threshold predictedby this quadratic interaction failure criterion is 33.4 ktaps, compared with 46 ktapsaccording to the experiments. Mode ifi threshold can be specified from the comparison ofthe tensile strain ratio and the shear stress ratio, which is shown in Fig. 4.22. As seenfrom the figure, after a large increase in the shear ratios, they remain quite constant past99C)0I-CuC,).QC) Co- 0CS 1CSCl)65431Ca010—1—265430 1 2 3 4 5Distance from the support On.)10 microsec. — 20 microsec. __ 30 microsec. 40 microsec.0—1—2—340 microsec.0 1 2 3 4Distance from the support On.)—_ 50 microsec. —_ 60 microsec.—a— 70 mlcrosec.5Fig. 4.20. Distributions of shear force along the half beam in the early stage of response.1003. 2.5C0E2—.— QICF (10 elem)_-_QICF (14 elem)_-_- LICF (10 elem)+ Mode I Exp. Results—A—- QICF (12 elem)A Shen&Jones’Fig. 4.21. Midspan permanent deformations vs. impulses for a 0.375 in. by 1 in. by 8 in. beamusing the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) andthe linear interaction failure criterion with shear stress from the FE analysis (LICF).0 10 20 30 40 50 60 70Impulse (ktaps)1011.10.6 -: she70Impulse (ktaps)__ tensile ratio (10)__- tensile ratio (12) —— tensile ratio (14)_- shear ratio (10) —. shear ratio (12) —h-- shear ratio (14)Fig. 4.22. Comparison of tensile strain ratios and shear stress ratios vs. impulses for a 0.375 in.by 1 in. by 8 in. beam using the quadratic interaction failure criterion wlih shear stress from theFE analysis (QICF).10261.6 ktaps. Then this is specified as the mode ifi threshold, which is very close to theexperimental result of 64 ktaps.The mesh is then further refined to 14 elements, with the first element in the 10 equalelement mesh subdivided into 4 equal length elements and the second element halved. Theresults in the mid-span permanent deformation are almost the same as from the 12 elementmesh. However, a large difference is observed from the comparison of tensile strain ratioand shear stress ratio (Fig. 4.22) because of the boundary layer effect (discussed inChapter 2). The distribution of the first order shear force at the midpoint of each elementalong the beam at 60 Lsec, which is the duration of the load, is shown in Fig. 4.23. As isevident from this figure, the shear force at the midpoint of the first element is smaller thanthat in the middle of the second element. This is because the midpoint of the first elementis inside the boundary layer, within which the shear force erroneously decreases from themaximum to zero at the very beginning boundary of the beam. To avoid this, an elementgreater than is preferred, to be conservative.Analyses have also been done for the case of a 0.25 in. x 1 in. x 8 in. beam subjected tovarious impulses, using the same interaction failure criteria. Results are shown in Figs.4.24 and 4.25. In this example, the mode II and ifi thresholds are predicted to be 21.8and 39.6 ktaps, which compares well with the 28.7 and 46 ktaps reported fromexperiments.3. Failure Criterion Using Theory ofPlastic Work DensityTo assess the accuracy of the failure criterion using sectional plastic work density inmodelling the failure of the blast loaded beam, a series of finite element failure analysesincorporating this failure condition have been done for a 0.375 in. x 1 in. x 8 in. beamsubjected to different impulsive loading conditions. The mid-span permanent103543—. 0 2a) .o C00°’h. 0a)a, I-iCl,0—20Fig. 4.23. Distribution of shear force at 60 microseconds along half a 0.375 in. by 1 in. by 8 in.beam subjected to an impulse of 17.8 ktaps for a forty element mesh.1 2 3 4 5Distance from the support (in.)1043• 2.5.-0220.50. QICF (10 elem)Shen & Jones’—..— LICF (10 elem)—--- QICF (12 elem)Fig. 4.24. Midspan permanent deformations vs. impulses for a 0.25 in. by 1 in. by 8 in. beamusing the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) andthe linear interaction failure criterion with shear stress from the FE analysis (LICF).30Impulse (ktaps)A Mode I exp. results1.10.90.80.7 -0.6 -0.5—0.4 -0.3 -0.2 -0.10_.— tensile strain ratio (10 elem)— shear stress ratio (10 elem)—.-- tensile strain ratio (12 elem)—a--- shear stress ratio (12 elem)Fig. 4.25. Comparison of tensile strain ratios and shear stress ratios vs. impulses for a 0.25 in.by 1 in. by 8 in. beam using quadratic interaction failure criterion wfth shear stress from FEanalysis.105tensile strain ratioshear stress ratio20 30 40Impulse (ktaps) 50 60 70106deformations obtained from the analyses are plotted versus impulses applied on the beam.Two different meshes have been employed to represent one half ofthe beam: 10 elementsof equal length, and 12 elements with eight of equal length close to the mid-span and fourof half size near the support. Results are presented in Fig. 4.26, along with theexperimental data. As seen from the graph, neither curve predicts the trend of mid-spanpermanent deformation well. Furthermore, there is quite a large difference between thecurves obtained from the two different meshes. Analysis based on the 10 element meshpredicts the mode II failure at 25.6 ktaps, while the 12 element mesh analysis predicts 21.8ktaps. This suggests that analysis using the failure condition of sectional plastic workdensity is mesh dependent. This is probably because the sectional plastic work density iscalculated at each spanwise Gauss point. The coarser the mesh, the fewer the Gausspoints, and thus there are fewer sections to be considered in the failure condition, and themore unlikely the beam will fail. To overcome this constraint, the failure criterion ischanged into comparing the plastic work done in the first element (which is checked to bealways the highest) to the critical value. The results from the analysis using the changedfailure condition employing a 10 element mesh are also shown in Fig. 4.26. The resultsare still unsatisfactory as the trend of the mid-span permanent deformation is notaccurately portrayed and the mode II threshold is predicted as early as 25.6 ktaps. Thisamendment did not improve the prediction probably because the stress state at theboundary of the beam undergoing large plastic deformation is three dimensional, and theone dimensional beam theory cannot accurately describe the beam behaviour near theboundary. The plastic work density calculated near the boundary cannot represent the realvalue exactly, making it difficult to use the failure condition based on plastic work densitytheory to predict the rupture of the beam.Although the failure criterion based on plastic work density does not predict the failureof the beam satisfactorily, a detailed investigation of the distribution of sectional plastic321.50.50107Fig. 4.26. Midspan permanent deformations vs. impulses for a 0.375 in. by 1 in, by 8 in. beamusing failure conditions based on sectional plastic work density (SPWD).. 2.5C0E0)Ca,Ca,Ea,0Ca,0U,0 10 20 30 40 50 60 70Impulse (ktaps).. SPWD (10 elem) SPWD (12 elem) A Experimental ResultsShen & Jones’ integrated SPWD108work along the beam will help in understanding the distribution of plasticity in the beamwhen subjected to different impulses. Shown in Fig. 4.27 are a series of spanwisedistribution of sectional plastic work density for a 0.3 75 in. x 1 in. x 8 in. beam subjectedto an impulse of 17.8 ktaps, which caused mode I failure in the experiments. The sectionalplastic work density shown is evaluated at the midpoint of each element. The samedistributions for the same beam subjected to 42.9 ktaps and 64 ktaps (which caused modesII and ifi failure experimentally) are presented in Figs. 4.28 and 4.29. In the figures, 107and 46 sec are the times of mode II failure predicted by the quadratic interaction failurecriterion for the two cases, respectively, while 120 and 66 J.Lsec are the times of mode IIfailure predicted by failure condition based on sectional plastic work density, respectively.From these figures it is obvious that the plastic work done in the first element is thehighest. The development with time of the plastic work in the halfbeam for mode I failurecan be clearly observed from Fig. 4.27. The plasticity first occurs close to the boundary,then develops about the middle of the half beam, and finally appears in the middle area ofthe (whole) beam. In mode II, rupture of the beam occurs before much plasticity developsin the middle of the beam (shown in Fig. 4.28). In mode ifi the beam fails even earlier andthe length of the beam where there exists significant plasticity is far shorter (refer to Fig.4.29). That is, as the load increases, the failure of the beam moves to mode II and thenifi, and the plasticity in the beam becomes more localized. This agrees with experimentalobservations.Another detailed study of the time histories at different energies gives some idea aboutthe different energetic characters of different failure modes. The time histories of externalwork, kinetic energy, and plastic work done in the half beam for a 0.375 in. x 1 in. x 8 in.beam subjected to three typical loading cases (17.8, 42.9, and 64 ktaps) which correspondto modes I, II, and ifi failure of the beam, are presented in Figs. 4.30 (a), (b) and (c)400100080070010015005001093000-x0:20000a.0C0a00)Distance from the iupport (en.)_. 25 microsec. ....÷.... 50 microsec.-__75 miorosec. 100 mioroseo.>.6000¶so300C0iS 200e0)Distance from the Lpport (in.)_.. 150 miorosec. ........ 175 microsec._-_200 microsec.toC10000Distance from the lupport (in.)225 microsec._,_250 microsec.___275 microseo. -__ 300 mlorosec.Fig. 4.27. Distributions of sectional plastic work density along half a 0.375 in. by 1 in. by 8 in.beam subjected to an impulse of 17.8 ktaps._-__325 microsec.2500>%20004-Cl)CCu01500C)CuQ-iooocuC0Cu5000—.--9 microsec.——75 microsec.—--25 microsec.—-- 100 microsec.—-_50 microsec._-— 107 microsec.110Fig. 4.28. Distributions of sectional plastic work density along half a 0.375 in. by 1 in. by 8 in.beam subjected to an impulse of 42.9 ktaps.Distance from the support (in.)0 1 2 3 41111500>‘0C1000Io 50013Cl)0_ 25 microsec.Fig. 4.29. Distributions of sectional plastic work density along half a 0.375 in. by 1 in. by 8 in.beam subjected to an impulse of 64 ktaps.Distance from the support (in.)6 microsec.__- 46 microsec.1591050 0.5 Time (millisk.)Fig. 4.30. Time histories of different energies for half a 0.375 in. by 1 in. by 8 in. beam subjectedto impulses of (a) 17.8 ktaps. (b) 42.9 ktaps. (C) 64 ktaps.11212001000600600Cw4002000a765[32020C.external workkinetic energyplastic work01.5 2113respectively. As seen from Fig. 4.30 (a), during the mode I failure the majority of theexternal work is absorbed into the plastic work in the beam; this corresponds to theplasticity fully developing in the whole beam as discussed above. In mode II (Fig. 4.30(b)) the beam still acquires quite an amount of kinetic energy after the plastic workabsorbs the majority of the external worlç this corresponds to the plasticity developed inthe beam decreasing. Finally, in mode ifi failure (Fig. 4.30 (c)) most of the external workis turned into kinetic energy of the beam. The plastic work in the beam is relatively verysmall; this corresponds to the beam failing at an early time and plasticity in the beam islocalized.The residual kinetic energy is further plotted vs. impulse for half a broken beam of 0.3 75in. by 1 in. by 8 in., as shown in Fig. 4.31. The highest velocity of the flying broken beamis 8636 in./s (220 m/s). Similar results are obtained for the 0.25 in. by 1 in. by 8 in. beam,which are shown in Fig. 4.32. In this case, the velocity of the flying broken beam can beas high as 13975 in./s (355 mIs), which is faster than the sound speed. Research on theblast loaded circular plates [29] also shows the same trend in the residual kinetic energyvs. impulse relation.4.4.3 Convergence and Energy ConservationHaving discussed the accuracy of the failure conditions in modelling the actual failure,the accuracy of the finite element calculation and the time integration is still of interest.This is to be discussed next.To check the effect of the time step on the accuracy of the time integration, three timesteps, 1, 0.5, and 0.1 sec, have been used in analyzing the beam of 0.375 in. x 1 in. x 81142015 -C.LIC0)5 1OCoW1E0C5-QIC0 I I I0 10 20 30 40 50 60 70Impulse (ktaps)QIC (10 elements) . LIC (10 elements)Fig. 4.31. Kinelic energy of half a broken beam of 0.375 in. by 1 in. by 8 in. subjected to differentimpulses.-C2>Cw0C115LIC3025 -20—•0J1510 -5-0—0 10 20 30 40 50 60 70Impulse (ktaps)QIC (10 elements) . LIC (10 elements)Fig. 4.32. Kinetic energy of half a broken beam of 0.25 in. by 1 in. by 8 in. subjected to differentimpulses.QIC116in. subjected to an impulse of 17.8 ktaps. A 10 equal element mesh is used. The mid-spanpermanent deformations obtained are listed as follows for comparison:Table 4.1 Convergence of mid-span permanent deformations with time step.Time step max Aw mm Aw Aw(.isec) m m m1 0.64495 0.566235 0.602350.5 0.64244 0.55977 0.601080.1 0.64034 0.55758 0.59896The difference between results from calculations using 1 and 0.1 sec as time steps isonly about 0.56%. The answer obtained using 1 j.sec can be considered accurate enough.To check the effect of element mesh on the accuracy of the finite element calculation,three element meshes have been used: 8, 10, and 16 equal elements representing the halfbeam, for the case of 0.375 in. x 1 in. x 8 in. beam applied with an impulse of 17.8 ktaps.The mid-span permanent deformations Aw are shown below:Table 4.2 Convergence ofmid-span permanent deformations with mesh.No. of elements max Aw mm Awin. m m8 0.63716 0.55488 0.5960210 0.64495 0.56235 0.6036516 0.64772 0.56528 0.6065From the table it can be observed that the results are converging well.117For the energy conservation check, two cases of a 0.375 in. x 1 in. x 8 in. beamsubjected to 17.8 ktaps and 52.9 ktaps are analyzed using different time steps. The resultsare presented in Figs. 4.33 and 4.34 respectively. In these figures the time histories of theenergy difference ratio is shown. Energy difference ratio rE is defined asEt -Ek plrE= E(4.4)extwhere Ee,, Ek, and E1 correspond to external energy, kinetic energy, and plastic workdone in the half beam. Since the Eext should be greater than the summation ofEk andwith the difference being strain energy of the beam, rE should always be greater than 0.However, when a larger time step is used, the numerical error can cause negative rEs. Asobserved from the figure, when rE is less than 0, it is always greater than -1.5%, which isstill acceptable. With the decrease in the time step, the rE is positive, so the resultsconverge to the correct value. Therefore, the energy is conserved when a small enoughtime step is used.4.5 SummaryIn this Chapter, after the series of experiments conducted by Menkes and Opat areintroduced, the numerical simulation of the experiments by using finite element analysisincorporated with the failure conditions are presented. The results obtained by employingdifferent failure criteria are compared to the data from experiments in the discussion. Theconvergence is also presented and, as shown, the results are converging as the time stepand mesh are refined.1180.10.09 -0.08 -o0.07 -4-0.06-C0.05 -a)0.04->0.03-C0.02 -0.01 -0—0.01 I I0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1Time (millisec.)—.-- ts=1 microsec. (10 elem)__ts=0.1 microsec. (12 elem)... ts=0.1 microsec. (10 elem)Fig. 4.33. Time histories of energy difference ratio using different time steps for a 0.375 in. by 1in. by 8 in. beam subjected to an impulse of 17.8 ktaps.119Fig. 4.34. Time histories of energy difference ratio using different time steps for a 0.375 in. by 1in. by 8 in. beam subjected to an impulse of 52.9 ktaps. The rupture of the beam is predicted bythe quadratic interaction criterion with shear stress from the FE analysis.0.020.015.2 0.01a)0.005a)Va)CW—0.005—0.01—0.0150 0.5 1 1.5 2ts=1 microsec. (10 elem)—i— ts=O.1 microsec. (10 elem)Time (millisec.)ts=0.5 microsec. (10 elem)ts=0.1 microsec. (12 elem)2.5120According to the comparison, the quadratic interaction failure criteria using shear stressfrom the finite element analysis provided the closest answer, but at the cost of using adenser mesh. However an extremely dense mesh should be avoided because of theboundary layer effect of the higher order beam theory solution. The quadratic interactionfailure condition using wall reaction force provides a good prediction as well whileallowing a comparatively coarse mesh to be used to obtain accurate results; thus thistechnique is more practical for engineering applications.From the discussion, the effect of interaction cannot be neglected, particularly in themode ifi failure. The membrane effect still plays an important role in the shear dominantfailure of the beam.Additionally, the ruptures of the beam predicted from the analysis are earlier than theexperimental results, which is desirable for engineering practice. This conservatism isprobably because the beam ends in the experiments may not be fi.illy restrained accordingto the configuration of the experiment and the length of plastic hinge cdi might be smallerthan in reality.121CHAPTER 5CONCLUSION5.1 SummaryA failure analysis procedure to predict the rupture of the impulsively loaded ductileclamped beams which exhibit geometric and material non-linearities has been presented inthis study.A detailed look into the finite element formulation of the problem is first presented.With the kinematics of the higher order beam theory and the non-linear strain displacementrelation by using the principle of virtual work and Total Lagrangian Approach, thegoverning variational equation of motion is obtained. It is yet to be noted that the beamkinematics is restricted to the range of small strains and moderately large displacements.After the finite element spanwise discretization is applied, the variational governingequation of motion reduces to a system of equations of motion which are differentialequations with respect to time. The constitutive model which the system of equations ofmotion are based on considers elastic-plastic isotropic strain hardening material behaviourby von Mises yield criterion and associate flow rule and strain rate sensitivity by CowperSymonds relationship. Then the differential equations are integrated over the time domainby the central difference method. Lumped mass matrix is used to improve the accuracyand efficiency. A conservative critical time step is determined for stability. The governing122differential equations of motion now further reduce to a system of algebraic equationswhich can be solved easily.Based on finite element method calculations interaction failure criteria and the failurecondition using plastic work density theory are proposed.The premise behind the interaction failure criteria is to consider the contribution of bothtensile tearing and transverse shearing effects to the failure of the beam. Thus it canexplain the overlap of mode II and ifi failures between the thresholds. The beam ismodelled as an elastic-plastic beam with a plastic hinge at the supports. The plastic hingeassumption is employed in calculating the tensile strain. The shear stress can be obtainedeither from wall reaction forces (from equilibrium) or first order shear force (from thefinite element analysis). Tensile strain ratio and shear stress ratio are used to calibrate theinfluence of the tearing and shearing effect, respectively. Two interaction functions havebeen investigated, namely quadratic interaction and linear interaction.The concept on which the failure criterion using the plastic work density is based is thata crack initiates at a point where the density of plastic work reaches a critical value. It isassumed that the beam does not fail until the whole section fails since the rupture alwaysoccurs in a section. Therefore the failure conditions can be put that the beam will not failuntil the sectional plastic work density reaches a critical value.Because when the rupture occurs there still remains quite a lot of residual energy, thepost failure analysis on a blown off free beam is presented. In this way, the mid-spanpermanent deformation from the analysis is comparable to the experimental results.123To assess the accuracy of the failure conditions, analysis using different filureconditions have been carried out to simulate the experiments conducted by Menices andOpat. Results obtained from these analyses have been compared to the experimental data.For the quadratic interaction filure conditions, encouraging results have been obtained. Itis to be noted that Menkes and Opat admitted in their paper that the specification of modeII and ifi thresholds are highly subjective. Therefore, the conservative results obtainedfrom the analysis are desirable in the engineering sense.As stated in the introduction, the objective of the present study is to determine a simplefailure model to be incorporated into the finite element analysis so that the onset of themode II and ifi failures of impulsively loaded beam can be predicted. This has beensatisfactorily accomplished. Furthermore, some insight into the characteristics of thedifferent failure modes is obtained from the discussion about the numerical results.5.2 Areas Requiring Further InvestigationIn the course of the present work, several assumptions are made without confirmation;therefore they might be potential sources of error. Some of these are listed as follows forfurther investigation.During the formulation of the finite element analysis of the problem, it is assumed thatthe strain is everywhere small. However, when the rupture of the beam occurs due toexcessive strain, the strain becomes significant. This might have some effect on thevalidity of simplification in the formulation.124In the present work, the rupture strain of the beam is chosen to be equal to the staticuniaxial rupture strain. The ultimate shear strength is determined based on the von Misesyield criterion and the ultimate tensile strength of the material. That the dynamic shearstrength appears to be much higher than the static one is taken into account by using a“conservative” number. It can be noted that the determination of these parameters isbased on several assumptions which are not confirmed. An investigation is thereforeneeded to check the sensitivity of the results to these two parameters.During the modelling of the impulsive loading condition, a rectangular pressure pulse isused and an assumption was made as to the duration of the load. Since the assumptionson which the load modelling is based are not well confirmed, this introduces the possibilitythat the modelling of the load shape and duration might have some influence on the failureof the beam. A parameter study is needed to confirm or improve the modelling of theblast loading condition in the future.While calculating the bending strain at the support, the length of the plastic hinge isdefined as l=ah, where cx, the plastic hinge parameter, is chosen to be two. However, asthe beam goes into the membrane stage, the length of the plastic hinge will be greater, thusdecreasing the bending strain and delaying the onset of the mode II failure of the beam,which is closer to the experimental data. With the increase in the impulsive loading, thelength of the plastic hinge will decrease, as indicated in the figure on distribution of plasticwork along the beam. This will have an effect on the time of rupture of the beam and thusthe mid-span permanent deformation. As discussed above, the specification of the lengthof the plastic hinge, hence the plastic hinge parameter a, will have an effect on theaccuracy of the interaction failure model; possibly refinements are necessary.125References1. Jones, N. “A Literature Review of the Dynamic Plastic Response of Structures,”Shock and Vibration Digest, Vol. 7, August 1975, pp.89-105.2. Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 1and 2, Shock and Vibration Digest, Vol. 10, September 1978, pp.21-33 and October1978, pp.13-9.3. Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 3,Shock and Vibration Digest, Vol. 13, October 1981, pp.3-16.4. Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 4,Shock and Vibration Digest, Vol. 17, February 1985, pp.35-47.5. Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 5,Shockand Vibration Digest, Vol. 21, 1989, pp.3-13.6. Menkes, S. B. and Opat, H. J., “Broken Beams”, Experimental Mechanics, Vol.13,1973, pp.480-486.7. Jones, N. “Plastic Failure of Ductile Beams Loaded Dynamically.”, Journal ofEngineeringfor Industry, Vol. 98, 1976, pp.13 1-136.8. Jones, N., “On the Dynamic Inelastic Failure of Beams”, Structural Failure, T.Wierzbicki and N. Jones, eds., Wiley, 1989, pp.133-59.9. FoIz, B. R., “Numerical Simulation of the Non-Linear Transient Response of SlenderBeams”, MA.Sc. thesis, University of British Columbia, Vancouver, Canada, April1986.12610. Fagnan, J. K., “Failure Prediction of Blast Loaded Ductile Beams UsingApproximate Deflection Profiles”, Ph.D. Thesis Proposal, University of BritishColumbia, Vancouver, Canada, June 1991.11. Shen, W. Q. and Jones, N., “A Failure Criterion For Beams Under ImpulsiveLoading”, International Journal of Impact Engineering, Vol. 12, No. 1, 1992,pp.101-121.12. Jones, N. and Shen, W. Q., “Criteria for the Inelastic Rupture of Ductile MetalBeams Subjected to Large Dynamic Loads”, Structural Crashworthiness andFailure, N. Jones and T. Wierzbicki, eds., Elsevier Applied Science, 1993, pp.95-130.13. Jiang, J., and Olson, M. D., FENTAB version 2, Department of Civil Engineering,University ofBritish Columbia, Vancouver, 1990.14. Washizu, K., Variational Methods in Elasticity & Plasticity, 3rd Edition, PergamonPress, 1982.15. Levinson, M., “A New Rectangular Beam Theory”, Journal ofSound and Vibration,Vol. 74, 1981, pp.81-87.16. Bickford, W. B., “A consistent Higher Order Beam Theory”, Developments inTheoretical andAppliedMechanics, Vol. 11, 1982, pp.137-50.17. Heyliger, P. R. and Reddy, 3. N., “A Higher Order Beam Finite Element for Bendingand Vibration Problems”, Journal of Sound and Vibration, Vol.126(2), 1988,pp.309-326.18. Reddy, J. N., “A simple Higher-Order Theory for Laminated Composite Plates”,Journal ofAppliedMechanics, Vol. 51, December 1984, pp.745-52.19. Reddy, J. N., Energy and Variational Methods in Applied Mechanics. New York:John Wiley. 1984.20. Palmer, A. C., StructuralMechanics, Clarendon Press, 1976.12721. Cook, R. D., Malkus, D. S. and Plesha, M. E., Concepts and Applications ofFiniteElement Analysis, 3rdF4ition, John Wiley & Sons, 1989.22. Karasudhi, P., Foundations ofSolidMechanics, Kiuwer Academic Publishers, 1991.23. Olson, M. D., CIVL 538 Course Notes: Advanced Topics in the Finite ElementMethod, University ofBritish Columbia, Vancouver, BC, January, 1993.24. Olson, M. D., CIVL 537 Course Notes: The Finite Element Method, University ofBritish Columbia, Vancouver, BC, September, 1992.25. Jiang, J., and Olson, M. D., Documentation for FENTAB Department of CivilEngineering, University ofBritish Columbia, Vancouver, 1990.26. Nonaka, T., “Some Interaction Effects in a Problem ofPlastic Beam Dynamics, Parts1-3”, Journal ofAppliedMechanics, Vol. 34, 1967, pp.623-4.27. Fagnan, J. R., FENTABF, Department of Civil Engineering, University of BritishColumbia, Vancouver, 1991.28. Humar, J. L., Dynamics ofStructures, Prentice Hall, 1990.29. Olson, M. D., Fagnan, J. R. and Nurick, G. N., “Analysis of the Deformation andTearing ofBlast Loaded Circular Plates”, to be published.