UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Failure analysis of blast loaded ductile beam Luo, Qing 1994

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1994-0530.pdf [ 2.82MB ]
Metadata
JSON: 831-1.0050411.json
JSON-LD: 831-1.0050411-ld.json
RDF/XML (Pretty): 831-1.0050411-rdf.xml
RDF/JSON: 831-1.0050411-rdf.json
Turtle: 831-1.0050411-turtle.txt
N-Triples: 831-1.0050411-rdf-ntriples.txt
Original Record: 831-1.0050411-source.json
Full Text
831-1.0050411-fulltext.txt
Citation
831-1.0050411.ris

Full Text

FAILURE ANALYSIS OF BLAST LOADED DUCTILE BEAMS  by Qmg Luo B.Sc., Dalian University of Technology, 1991  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Civil Engineering)  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA September 1994 © QingLuo, 1994  _  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  (Signature)  Department of  / V/  L  &A;/ 7/ A1  The University of British Columbia Vancouver, Canada Date  DE-6 (2/88)  CA  -  /  ‘  //IY  II  Abstract  A failure analysis procedure based on the elastic-plastic finite element analysis of the higher order beam is presented to predict the mode II and ifi failure of impulsively loaded ductile clamped beams.  The variational equation of motion of the problem is first obtained by using the principle of virtual work and the Total Lagrangian Approach together with the kinematics of the higher order beam theory and the nonlinear strain-displacement relation.  The Finite  Element Method is employed to discretize the beam spatially initiating the numerical solution procedure.  The constitutive model considers elastic-plastic isotropic strain  hardening by von Mises yield criterion and associated flow rule and strain rate sensitivity by Cowper-Symonds relationship.  Equations of motion are integrated by central  difference method in the time domain.  Based on the finite element simulation, two failure criteria are proposed.  In the  interaction failure criteria, contributions from tensile tearing and transverse shearing to the damage of the beam are considered by including tensile strain ratio and shear stress ratio in the criteria. Both linear and quadratic models are investigated. The beam is modelled as an elastic-plastic beam with a plastic hinge at each support to calculate the total strain at the support.  The shear stress comes either from the wall reaction obtained from  equilibrium or directly from the finite element analysis of the higher order beam theory. In the sectional plastic work density criterion, the beam is assumed not to fail until the sectional plastic work density exceeds the critical value.  Post failure analysis is also  included to take account of the residual energy of the beam at failure.  111  Numerical simulations have been carried out for the experiments of blast loaded beams. Comparison with the experimental results suggests the quadratic interaction criteria is suitable for this type of analysis.  iv  Table of Contents  Abstract  ii  Table of Contents  iv  List of Tables  vi  List ofFigures  vii  Acknowledgments  xi  1.  2.  3.  Introduction  1  1.1. Background  1  1.2. Purpose and Scope  6  Finite Element Formulation and Solution Procedure of the Problem  8  2.1. Introduction  8  2.2. Kinematics of the Higher Order Beam Theory  9  2.3. Variational Equation ofMotion  15  2.3.1. Variational Equation ofMotion  15  2.3.2. Governing Differential Equation and Boundary Conditions  19  2.4. Finite Element Formulation of the Equations of Motion  21  2.5. Modelling the Material Behavior of the Beam  29  2.6. Solution Procedure Over the Time Domain  39  2.7. Summary  43  Failure Analysis of the Problem  45  3.1. Introduction  45  3.2. Interaction Failure Criterion  46  V  4.  3.2.1. Calculation of the Influence of Tensile Tearing  46  3.2.2. Measure of the Influence of Transverse Shearing  54  3.2.3. Two Interaction Failure Criteria  58  3.2.4. Another Two Interaction Failure Criteria  59  3.3. Failure Criterion Using Plastic Work Density  60  3.4. Post Failure Analysis  64  3.5. Summary  65  Numerical Examples and Discussions  66  4.1. Introduction  66  4.2. Menkes and Opat’s Experiments  67  4.3. Modelling the Problem  69  4.4. Numerical Results and Discussion  72  4.4.1. Numerical Simulation of the Response Process  72  4.4.2. Discussion on Failure Conditions  80  4.4.3. Convergence and Energy Conservation  5.  113  4.5 Summary  117  Conclusions  121  5.1. Summary  121  5.2. Areas Requiring Further Investigation  123  References  125  vi  List of Tables  4.1  Convergence of mid-span permanent deformations with time step.  116  4.2  Convergence of mid-span permanent deformations with mesh.  116  vu  List of Figures  1.1  Three failure modes associated with the clamped metal beams loaded impulsively.  2.1  3  Definitions of positive transverse displacement of the neutral surface and rotation of a cross section of the beam at the neutral surface.  2.2  13  Shear force edge effect for a cantilever beam loaded at x=0 with a tip load P.  22  2.3  An idealized bilinear model of the elastic-plastic material.  31  2.4  Idealization of elastic-plastic linear strain hardening material behavior with strain-rate sensitivity.  33  3.1  Typical deflection profile as determined by FENTAB.  48  3.2  Plastic hinge in a thuly clamped beam.  50  3.3  Shear distribution along the length of the half beam obtained fromFENTABv.2.  55  4.1  Schematic drawing of experimental configuration.  68  4.2  Experiment results for series of 6061-T6 beams (0.25x1x8-in. beam).  70  4.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.  4.4  Time history of central deformation for a 0.375 in. by i in. by 8 in. beam subjected to an impulse of 52.9 ktaps.  4.5  73  75  Comparison of time histories of mid-span deformation obtained from the FE analyses including shear (FENTAB v.2) and not including shear (FENTAB).  4.6  Comparison of time histories of mid-span deformation during  77  vm  early stage of response obtained from the FE analyses including shear (FENTAB v.2) and not including shear (FENTAB). 4.7  Comparison of beam proffles during the early stage of response obtained from FENTAB and FENTAB v.2.  4.8  78  79  Mid-span permanent deformations vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion (QIC) and the linear interaction failure criterion (LIC).  4.9  81  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 using the quadratic interaction failure criterion (QIC).  83  4.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 using the linear interaction criterion.  85  4.11 Comparison of pull-ins vs. impulses obtained from a 0.375 in. by 1 in by 8 in. beam using the quadratic interaction criterion (QIC)  and linear interaction criterion (LIC) with a ten element mesh.  86  4.12 Comparison of times of failure vs. impulses for a 0.375 in. by 1 in by 8 in. beam using the quadratic interaction criterion (QIC) and linear interaction criterion (LIC) with a ten element mesh.  87  4.13 Mid-span permanent deformations vs. impulses obtained from a 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction criterion and the linear interaction criterion with a ten element mesh.  89  4.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 using the linear interaction criterion. 4.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 using  90  ix  the quadratic interaction criterion.  91  4.16 Time histories of total strain at the support for different finite element meshes. 4.17 Time histories of wall reaction for different finite element meshes.  92 93  4.18 Comparison of time histories of F and wall reaction for an eight element mesh.  96  4.19 Comparison of time histories of F and wall reaction for an twenty element mesh.  97  4.20 Distributions of shear force along the half beam in the early stage of response.  99  4.21 Mid-span permanent deformations vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) and the linear interaction criterion with shear stress from the FE analysis (LICF).  100  4.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 quadratic interaction failure criterion with the shear stress from the FE analysis (QICF).  101  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.  103  4.24 Mid-span permanent deformations vs. impulses for a 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) and the linear interaction failure criterion with shear stress from the FE analysis (LICF). 4.25 Comparison of tensile strain ratios and shear stress ratios vs.  104  x  impulses for a 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion with the shear stress from the FE analysis (QICF).  105  4.26 Mid-span permanent deformations vs. impulses for a 0.375 in. by 1 in. by 8 in. beam using failure conditions based on  sectional plastic work density (SPWD).  107  4.27 Distribution of sectional plastic work density along half a 0.3 75 in. by 1 in. by 8 in. beam subjected to an impulse of 17.8 ktaps.  109  4.28 Distribution of sectional plastic work density along half a 0.3 75 in. by 1 in. by 8 in. beam subjected to an impulse of 42.9 ktaps.  110  4.29 Distribution 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.  111  4.30 Time histories of different energies for haIfa 0.375 in. by 1 in. by 8 in. beam subjected to impulses of(a) 17.8 ktaps (b) 42.9 ktaps (c)64ktaps.  112  4.31 Kinetic energy for haifa broken beam of 0.375 in. by 1 in. by 8 in. subjected to different impulses.  114  4.32 Kinetic energy for haifa broken beam of 0.25 in. by 1 in. by 8 in. subjected to different impulses.  115  4.33 Time histories of energy difference ratio using different time steps for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 17.8 ktaps.  118  4.34 Time histories of energy difference ratio using different time steps for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 52.9 ktaps.  119  xi  Acknowledgments  This work is possible due to the efforts of many people. Particularly, Professor Olson has given me advice, supervision, and guidance which have been invaluable. Also, I would  like to thank Dr. J. Jiang, Mr. Nagaraja Rudrapatna and Mr. Julien Fagnan for their helpfi.il discussions. Special thanks to Mr. Scott Tomlinson for his patient typing, proof reading,  and moral support.  Financial assistances were provided by the graduate research  assistantship from the Natural Sciences and Engineering Research Council grant and the University of British Columbia Graduate Fellowship.  Finally, I am grateful to my parents for encouraging me to study in Canada. Without them none of this would have been possible.  1  Chapter 1 Introduction  1.1 Background  It is becoming increasingly necessary to examine the influence of impact loading on structures in the design of various engineering systems. An understanding of the inelastic response of structures subjected to dynamic loads which cause permanent displacements or structural damage is required for various branches of engineering. For example, to avoid the destructive action of earthquakes on buildings, such an understanding can be used to guide the development of rational design procedures. This knowledge is also essential 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 vehicles subjected to blast.  Examining the influence of impact loading conditions on complete structures is quite complicated because the influence will depend on the type of structure on which the impact load is applied.  As a simplification, research is restricted to estimating the  influence 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 concerned with the inelastic response of beams subjected to impulsive loadings which cause permanent deformation or damage.  2 As 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 the impact loads are severe enough this transverse propagation of stress waves can cause failure of the beam in the form of spalling. This type of failure occurs very early (usually within microseconds of the initial impact) while the overall structural failure lasts longer by a few orders of magnitude.  It is therefore customary to uncouple this early wave  propagation behaviour of the beam from the overall structural response.  The present  work focuses on the overall response of the beam by assuming the external dynamic load is applied to the mid-surface of the beam instantaneously.  Menkes and Opat [6] conducted a series of experiments on the overall inelastic structural response of filly clamped metal beams with rectangular cross sections which are subjected to transverse impulsive loadings. From the experiments, three failure modes were identified as the load is monotonically increased, shown in Fig. 1.1. The failures were 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 shear failure 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 and ifi 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 the  influence of finite displacements.  This analysis, which accounts for the influence of  membrane and bending forces, gave good agreement with the Menkes and Opat’s experimental data from the beams which suffered mode I failure. To predict the onset of  3  (a)  (b)  (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 at supports. (c) mode ifi, transverse shear failure at supports.  4 mode 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 threshold condition for a mode II failure can be obtained by equating the maximum strain to the static uniaxial rupture strain of the material.  A square yield criterion relating the  membrane force and the bending moment is used. To estimate the threshold for mode ifi failure, the concept of transverse shear slide (which is analogous to the concept of a plastic bending hinge) is specified as an idealization of rapid changes of slopes across a short length of the beam. The mode ifi failure is then assumed to occur when the amount of the transverse shear sliding at the supports reaches the beam thickness. Here, a square yield condition relating the transverse shear force and bending moment is used. It is stated that encouraging agreement with the experimental data for the onset of both mode II and ifi has been obtained. However, Jones neglected elastic deformation and the effect of transverse shear on deformation. The interaction effect between the tensile and shear action (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 when the external dynamic energy is significantly greater than the maximum amount of strain energy 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 energy must be at least five times larger than the corresponding strain energy which can be absorbed in a wholly elastic manner. Moreover, when the durations of external pressure pulses are comparable to the natural period of elastic vibration, the material elasticity can also have a significant influence on the response of the beam.  Fagnan [10] predicted the onset of mode II and III failures while accounting for the effect of finite displacement based on the finite element procedure for the classical Euler Bernoulli beam proposed by Folz [9]. To predict the mode II failure, Fagnan also used the  5  plastic hinge model for calculating the maximum strain at the support, and defines that mode II failure occurs when the maximum strain reaches the rupture strain in the static uniaxial tensile test. For the mode ifi failure, it is assumed that rupture of the beam occurs when the shear stress exceeds the ultimate shear strength. Because analyses based on the classical Euler-B ernouffi beam cannot estimate shear stress, it is calculated from the wall reaction force obtained from the equilibrium of the symmetric beam. Yet, this theory still does not consider the transverse shear effect on the deformation of the beam and the interaction between the tensile and shear action. Furthermore, the Euler-Bernoulli beam theory is proved to be fundamentally inadequate for the analysis of impact conditions. Nonphysical results have been reported that waves of infinitesimal wavelength (which means discontinuities) propagate with infinite velocity. It can be shown that satisfactory results are achieved only for lower modes and the shear-stress series for this theory are non converging [6].  According to Ref. [1], the transverse shear effect is believed to have a more important influence 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 an important role in theories of beam vibration and dynamic behaviour under impulsive loading and it should not be ignored in the analysis. As to the interaction, no sharp distinction is found between failure modes II and ifi. Failures involving both of these modes were observed when the beams were subjected to impulsive loadings in the range between mode II and ifi thresholds. Therefore, interaction effects of both tensile and shear actions should be included in the failure conditions.  A more universal energy criterion to predict the inelastic failure modes of beams loaded impulsively and which includes effects of bending, membrane and shear action, is suggested by N. Jones and W.Q. Shen [11,12]. The rupture of the beam is assumed to  6 occur 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 work to 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, this method assumes a rigid-plastic material, which neglects the elastic deformation.  An  empirical relationship between the plastic hinge length parameter and the energy ratio 13 is not confirmed. Further experimental and theoretical studies are required to assess the accuracy of the failure condition.  The present research overcomes many of the above limitations. As explained later, the present work is based on elastic-plastic finite element analysis [13]. It considers the non uniform transverse shear effect by using a higher order beam theory. And, it includes the interaction effect of tensile and shear actions in the failure conditions.  1.2 Purpose and Scope  As stated in the previous section, the purpose of the present work is to develop a simple failure criterion to predict the onsets of mode II and ifi failure for impulsively loaded beams 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 be solved.  Only the overall structural response of the beam is considered, which includes large inelastic deformation and/or rupture.  Within this research scope, initially straight  7 rectangular beams, with at least one axis of cross-sectional symmetry, undergo moderately large planar deflection but small strains.  The finite element formulation using the higher order beam theoiy is presented in Chapter 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 stress state is two dimensional, von Mises yield condition and associated flow rule as well as Hooke’s law is used in the elastic-plastic constitutive relation. The Virtual Work Principle and 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. The spatially discretized equations of motion are then solved by integration over the time domain with central difference. Chapter 3 discusses the modeling of the failure prediction of the problem. Two kinds of failure criteria are formulated. In the interaction failure conditions, interaction is exhibited in the way that the failure function is a mathematical function of both tensile and shear contribution. In the failure criterion using the plastic work density theory, interaction is considered in that the sectional plastic work density includes contributions from both tensile and shear stresses. To evaluate the accuracy and efficiency of the failure modeling, numerical simulations of the experiments conducted by Menkes and Opat are carried out. Results obtained from the simulations are compared to the experimental results, and these are presented in Chapter 4.  Finally, Chapter 5  summarizes the research work, draws conclusions from the comparisons, and outlines areas requiring further investigation.  8  Chapter 2 Finite Element Formulation and Solution Procedure of the Problem  2.1. Introduction  In order to study the failure of the blast loaded beam, it is necessary to know the displacement field, stress, and strain state of the beam. This is accomplished by using the Finite Element Method.  The solution algorithm, which is capable of numerically  simulating the transient dynamic response of the higher order beam undergoing geometric and material non-linear behaviour, is presented in this chapter. First, three different beam theories 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 equations of motion. After applying the Finite Element approximation of the displacement field, the Finite Element formulation of the equation of motion is obtained.  The elastic-plastic  material model is also employed, taking account of plasticity by von Mises yield criterion and associated flow rule. Finally, the transient response of the beam is arrived at by integrating the equations of motion using the central difference method.  Although a similar formulation can be found in Ref. [9], the higher order beam instead of the Euler-Bemouffi beam is employed in the present work to specify the beam kinematics. 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-  9 strain relation in the present work. Moreover, the external damping is now taken into account in the equation of motion.  2.2 Kinematics of Higher Order Beam Theory  The following solution procedure is only relevant to the bending effects of a slender beam undergoing large deformations.  Large deformations are defined having the  maximum lateral deflection experienced by the beam large relative to its height, although the deflection is small relative to the longitudinal dimension of the beam. The torsion-free bending of the beam can be realized in the (x-z) plane by placing certain geometric and loading restrictions i.e. the cross-section of the beam has (x-z) plane as its longitudinal plane of symmetry and the resultant of transversely applied loads lies in this longitudinal plane of symmetry. For a rectangular beam of interest here, the x axis is set to coincide with the centroidal axis of the beam and the z axis is chosen to be the transverse symmetric axis of the beam cross section. Since the longitudinal dimension of a slender beam is much larger than its lateral dimensions, it is assumed that the stress components a 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 theory assumes that cross sections normal to the undeformed neutral axis remain plane, undistorted, and normal to the deformed neutral axis. The shear deformation is therefore neglected. Following this assumption, as well as those previously made, the displacement field can be obtained  10  (x,y,z)= u(x)— z 1 u  ôw(x) (2.2)  (u 2 x,y,z)=0 3 (x,y,z) = w(x) u  where u(x) and w(x) represent, respectively, the horizontal and transverse mid-surface displacements of the beam. Shear deformation is accounted for in the Timoshenko beam theory [19,20], which assumes planes which are normal to the beam axis in the undeformed state remain plane in the deformed state but no longer normal to the beam’s deformed neutral axis. The displacement field obtained is  u (x,y,z) = u(x) + z(x) 1 (x,y,z)=0 2 u (x,y,z)= w(x) 3 u  (2.3)  where (x) is the rotation of the cross sectional plane. However, the relaxation of the classical Euler-Bernoulli assumption in Timoshenko beam theory leads to the contradiction that the resulting stress field no longer satisfies the shear-free boundary condition on the surfaces of the beam.  A correction factor, the Timoshenko shear  coefficient, is necessary to correct the contradictory shear stress distribution over the cross section 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 parabolic distribution of the transverse shear strain, so there is no need to use the shear correction coefficients. The higher order beam theory correctly accounts for the stress free boundary conditions on the upper and lower surfaces of the beam.  As observed from the Timoshenko beam theory, the longitudinal displacement is only a linear expansion with respect to the thickness coordinate. In the current higher order beam theory, the displacement field is chosen to be of a more general form including higher order terms. Refs. [15] and [18] give a detailed description about this theory.  11 Because of the condition that the transverse shear stresses vanish on the upper and lower surfaces of the beam and be non-zero elsewhere, the use of a displacement field in which the longitudinal displacement is expanded at least as cubic functions of the thickness coordinates is required. Since the beam is slender, both the Poisson’s ratio effects and stress components other than longitudinal normal stress and transverse shear stress are neglected. The transverse deflection of the displacement field can be set constant through the beam thickness. That is, the displacement field begins with,  (u 1 x,y,z) = u(x) + zNJ(x)+ z (x) + z 2 4(x) 3 (x,y,z)=O 2 u (x,y,z)= w(x) 3 u where  v is the rotation of a normal to the axis of the beam.  (2.4)  The functions (x) and 4(x)  are determined from the shear free condition on the top and bottom surfaces of the beam i.e. t(x4)=O  (2.5)  This is equivalent to requiring the corresponding shear strain to vanish on these surfaces, giving 1  &3  8W  3 z 2 8=-j---+----=N!+ z + 2 P+--  Setting  (2.6)  (x4) = o gives  (x)=O and  x)=—---(w+)  (2.7)  Therefore, the displacement field for the higher order beam theory becomes, if time dependence is taken into account,  12  2 4z  (x,y,z) = u(x,t) + z[W(x,t) — () 1 u  (w+  8w  (2.8)  (x,y,z)=O 2 u 3 (x,y,z) = w(x,t) u  This is with the positive transverse displacement of the neutral surface and the positive  rotation of a cross-section of the beam at the neutral surface, shown in Fig. 2.1, and the definition of positive longitudinal displacement same as for the Euler-Bernoulli beam theory. Indeed, 4(x), which is  (2.9) can be called the warping function since it reflects the extent of the deviation of the deformed beam cross section from an original plane surface.  The displacement field  implies that, in the higher order theory, the cross section is not only allowed to rotate relative to the neutral axis (as in the Timoshenko beam theory) but also to warp into a non-planar section, which is specified to satisfy the shear free boundary conditions at the top and bottom of the beam.  Since the beam undergoes finite deformation, Green’s strain tensor is employed to describe the state of deformation of the beam. With the displacement field established, the associated Green strain tensor (with reference to a Cartesian coordinate system) can be obtained by applying 1  ãuj  ôu.  ôu ôu  (i,j= 12,3)  (2.10)  to the higher order beam kinematics. Notice the restrictions imposed by assumptions discussed previously. Non-zero Green strain tensor components follow as  13  ‘(x,t) + -dx  w(x,t)  0  Fig. 2.1. Definitions of positive transverse displacement of the neutral surface and rotation of a cross section of the beam at the neutral surface. From Ref. [15]  14  +!()2  =  (2.11) £13 =(N’+)[1—4() ] 2  Eq. (2.11) can also be put into this form:  2 4z(z’  6ii=e+zKi+--)  2  (2.12) 83  1 =[1_4()2]q  ÔU  where the axial stram e = - + =  2 laW  (-i-) the first curvature id ,  —(-k + —j-). and the effective shear stram ( = i41+  =  -, the second curvature  6w -.  This is the non-linear strain-displacement relation for the higher order beam theory. As  observed from the relation, the von Karman strain is included to describe the geometric non-linearity caused by moderately large deflections of the beam.  The effect of von  Karman strain is realized to be important when the transverse deflection is equal to or larger 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 is everywhere small, i.e. (2.13)  15 2.3 Variational Equation of Motion 2.3.1 Variational Equation of Motion  With the higher order beam kinematics specified previously, this section will establish the variational equation of motion. Since the beam undergoes finite deformation and the Green strain tensor is still small everywhere in the beam, a Total Lagrangian Approach is adopted. In this approach, the rectangular Cartesian coordinates are fixed in space as the references of the body before and after deformation, and the coordinates defining points of the original undeformed body are employed for locating the points of the subsequently deformed body. For this problem, the x axis of the fixed rectangular Cartesian coordinates is set to coincide with the centroidal axis of the beam and the z axis is chosen to be the transverse symmetric axis of the beam cross section in the undeformed configuration, as discussed in the previous section.  The principle of virtual work is an avenue which leads to equations that govern the dynamic response of the structure.  For the present dynamic problem with finite  deformation 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 above expression,  °p,° lcd  ,°  V,° S represent, respectively, the mass density, the material  damping parameter, the volume, and the surface area (on which the surface traction is specified) of the beam before deformation. Sjj is the second Piola-Kirchoff stress tensor, u is the displacement vector, T is the prescribed surface traction vector, 6u is the small  16 virtual displacement vector which satisfies both compatibility and essential boundary conditions, and  the corresponding Green strain tensor. Each of the above terms is  specffied with reference to the undeformed configuration so as to be consistent with the Total Lagrangian Approach concept.  Here, the self weight of the beam has been  neglected 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 in the virtual work expression (2.14) is  .1 Su6EudV  =  ov  13 )dV & 13 11 + 2S & 11 .1(S  (2.15)  ov  because 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 Eulerian stress tensor Gjj (in which stress resultants are expressed) by the transformation equation  °p S, =—[a  ôu•  ôuãu.  (2.16)  (i,j,k= 123)  —  where p is the mass density in the current deformed configuration, and ô is the Kronnecker delta. For the present problem uj (i1,2,3) is set as for the higher order beam theory discussed in the previous section (Eq. 2.8). Because the Green strain tensor is everywhere small i.e. °p  <<1,  ---  (1,k=1,2,3) is believed to be far less than 1 and  11 p is considered to be valid. Therefore, S 11 i-’& )dV 13 I SejdV= f(a& ov  Introducing stress resultants,  ov  =  a, S 13  3 a  (2.17)  =  t,  i.e.  17  71= SodA ?Th  =fZ  V= ft[1_4()2]dA  m;=  S  ()2d4  j2 (2.18)  where 71 is the axial force, ml is the first order moment about the neutral axis, 771) is the third order moment about the neutral axis, V is the effective shear force, t is the first order shear force, and t is the second order shear force. Obviously, V  =  V  +  t. Here, °A  represents the original cross sectional area of the beam. Inserting the Green strain tensor components and stress resultants into Eq. (2.17) gives  +  =  fl-ö(-) + m 8() 1  +  —  V(NJ+ .)jix  +  (2.19) Eq. (2.19) can also be written as follows, using Eq. (2.12)  Ssjoe,  =  oic + m f[noe + 1  +  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), after substituting the higher order beam kinematics (Eq. 2.8), is as follows [16,17]  18  J °püöudV Dy  =  u 3 f(°puiöui+°pu  )dV  Ov  168i’  =Jf p Auöu+j p INJ&V—j  p  I--NI—j  p  8w 18’ 8w Iö(--)+pAw]dx Iwô(--)-Fj p  °L  (2.21) where  01  =  A similar  —bh, the moment of inertia 1 -  of  12  expression  the original beam  cross  section.  can be written for the third damping term in Eq. (2.14).  The work done by the external forces is  I TöudS  =  (2.22)  f[fou + fow]dx  Os  where f and f are the distributed axial and  transverse  loads acting on the beam with units  of force per unit length.  Finally, summarizing the above discussion,  o() 1 J[no(.) + flö() + m +°  p°Aüöu  +°  p°A5w  68  ..  p°Iöi.j,  —  +  +  Vö(j, +  8i 16 1 0’ 8w 8w 16 —j° p°I --ö —-j-° p°Iö(--) +-j° p°I--ö(--) ..  16 16 8w 8w 68 1 °Aüöu j°1’d °IVöWjj°1Cd IöWj1cd °INJö() d 8*8w 1 d °A*6w d °Ijö(j)-1-° 1 1 j .  —f8u  —  fow]dx =0  .  (2.23)  This is the variational equation of motion for the higher order beam theory. Starting from this equation, a complete set of governing differential equations of motion, natural  19  boundary conditions, and essential boundary conditions can all be derived for the new  beam theory, which is briefly discussed next. However, more important is that from now on the displacement field can be approximated by the finite element version and the algebraic equation system for the unknown variables can be derived to solve the problem numerically.  This is significant as the analytical solution is almost impossible for the  present problem, involving both geometric and material non-linearities.  2.3.2 Governing Differential Equations and Boundary Conditions  With the variational form of the equation of motion established, the governing differential equations will be persued in this section. This will be helpful when applying boundary 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 the coefficients of &i, ow, and 5j,, the governing differential equations are:  PAU+ICdAU  =  3 am —--  o  (71  —  O7l3  Ow -  —i- + + V) +  V  68 =  ..  16  ãi  68  .  16  (—iöPIw+ jI--) + (—-jKdIw+ j’d’)  f = pAM’ +  16  p1  Ofr -  1 -  O  16  0.j,  p1 —i- + )CdAW + j d’ E  1 -  (2.24) The appropriate boundary conditions are as follows:  20  natural boundary conditions  essential boundary conditions  71  u ‘I’  16 8w 8?712 1 1 8w 16 8W fl————+V--—pI+—pI——-—1cI\+—1cjI— 105 21 ax 21 105 ax .  W  8w ãx  Stress resultants are defined as in Eq. (2.18).  As a consequence of generalizing the displacement field beyond the traditional assumption (i.e. plane sections remain plane), the force resultants in the current higher order theory are also more generalized. Third order bending moments and second order shear forces occur which do not exist in the classical beam theory [161.  According to the results of variational principle, since the force boundary condition is unknown at the clamped boundary, u  w=  =  combinations of boundary conditions to verify this.  N’= 0.  Bickford [16] did three  The same three combinations of  boundary conditions were checked in the present work: (1) u = w =  u=w=V=0, u= w= w,  8w  =  and  and  (3)  u=w===0.  Both  results  =0, (2)  confirm  that  =0 is appropnate for a clamped end. The reason for this is that since u, .  w on the centerline  of the beam are used to descnbe the deformation of the  whole beam with the higher order beam kinematics (2.8), specifying uiO across the beam depth at the clamped boundary makes it equivalent to forcing u = w = satisfy ui=0 at every point on the clamped section.  =  =0 to  21 For a simple supported boundary, from the variational principle, the essential boundary condition should be uw=O. This is also confirmed by test calculations with combinations of different boundary conditions.  Also from the variational principle, the wall reaction force at a clamped boundary for a static beam will be  8w  3 am  (2.25)  Rw=fljj+V  Actually this is also the transverse reaction force at the end of a static beam of any boundary condition. If the von Karman effect is neglected  3 am  F=—---+V  (2.26)  is the appropriate transverse shear reaction force transmitted along the length of the beam for this theory. This is verified by Bickford [16]. Furthermore, he found the boundary layer 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 at the clamped boundary (x=L) V) approaches to 0 quickly, the higher order resultant  3 m  increases rapidly to supply the total resistance, and the transverse shear force transmitted along 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 Motion  Having obtained the variational equation of motion, to solve the present problem  involving both material and geometric non-linearities the Finite Element Method will be used to approximate the solution numerically. As usual the finite element approximation of the displacement field discretize the spatial domain of the beam into a number of  22  F  P  IL  Fig. 2.2. Shear force edge effect for a cantilever beam loaded at x0 with a tip load P. From Ref.[161.  23 subdomains, or finite elements. Within each element the approximate displacement field is assumed in terms of generalized nodal coordinates. According to the variational statement in Eq. (2.23), the transverse displacement should be at least twice differentiable and C  continuous, and the axial displacement u and rotation  w  should be at least once  differentiable and C° continuous [17]. Therefore, for a typical element, the displacement field u is approximated by: (2.27) 2  2  u(q,t) =  4  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 , 11 A =w and 3 A 2  for a typical element. Nj are the linear Lagrangian interpolation  functions and Mj are the Hermitian cubic interpolation functions.  (ri)=f 2 N  Ni(rU=1—f  =  =  131:11) +2(f)  (J  —  2(f)  (TI)= 1i[1_2(f)+(f)] 2 M  (2.29)  M4(fl) =  where 4 is the length of the element. As can be seen, the displacement field u is a function of both space and time, the shape functions are functions of space only, and nodal variables are functions of time only.  In this way the local separation of variables is  accomplished. In vectorial form, the appropriate set of generalized nodal coordinates for each element is:  1 deUi,NJiWi()  1 X  iT  ow (J , ,u 2 w 2 1V 1 2  ).  (2.30)  oui  )  1e ) iv e 1  ‘e I1re  1e J1 1ive ue  .  Le , Q  —‘  we  N  e 1 iir  1= ttl Ne  —  £  J  pU  {o’o’ o’11-’o}= n {io’o’ ‘io’o’-}= apt?  (t’€z)  z  =  i  dJ_  a1a. +  =  Jj  T[Z)  (EEc) Zi  (j  :si (j ) uO!Wnb uI14s uJ9 ) U9W9j joidA g joj pjj wowodsip UOTSJA 1UUI 1TU Otfl &irrnnsqn  soup1oo pozqiou jpou pajeTooss iosu uiis ua  tp  Jo suu u pssaidxo q u  UUT tp  ‘potsgqis ou pjj uwJds!p puTnss  xuww uowunj odqs aP[N]  (c)  =  s! [N]  iAE  fl  :sATwAuop UIT1 OM1 1SITJ S1J  ap[N]....ap  (wz)  :swoaq U!  uwp jdAi  o 0  n o 0  iui.up u  qii  11 Zp  0 tN  N 1 0  o  o  0 0  ‘iv 0 0  0 TN  T N 0  o  o  ioj pjj uwo1Jds!p  munb  oss  UAT  pu  14t1 —Ih--fl  L’J  utjj, qsow  a(.)  uondusqns  tj  25 00 00 0  0  0 0 B-  0  0 0 1 8M aM 1  0 0 1 aM 8M 2  iian 2 ôM 6M 1  ô’in 2 aM M 2  a1afl  a 1 ai  0  0  0 3 1 ôM 8M  0 3 2 ôM ôM  aa  ao  4 8M 8M 1  4 8M ôM 2  aaq  ao  0  00 0  0  0 0  0  00 00  0 1 ÔM  0 0 1 ôM M  a aM aM 2 3  thiô 4 aM ôM 2  0 0  0  0  00  0 3 aM ôM 3 aa 3 4 ôM 8M  0 M 8M 3 4 aa 4 ÔM aM 4  0  0  0 0  0 0 0 0  aa  Since the beam is discretized into elements and the displacement field satisfies the necessary continuity, the virtual work statement (2.14) for the discretized beams turns into: e=1  f  (2.36)  fTdudS  (s,Jö6,+°pööu+icdüöu)dv=  e=1  where NE represents the total number of elements in the beam.  For each element the work done by the internal forces:  =  where Ve aiid  6 are,  J[n& +  +  2 + V&p]dii moic  (2.37)  respectively, the volume and the length of the element. With e,  1C2, and 4i defined previously in Eq. (2.34):  I SjEj  =  6deTI[flBi + ?mB 3 2 + mB  +  VB  +  fSU&,ödeTrernt  flBde]d (2.39)  Ve  where reul#  =  1 J[nB  +  2 B 1 m  +  3 mB  +  VB  +  flBde]dll  (2.40)  (2.38)  , 1 K  26 is the internal force vector for a typical element which takes into account the internal  resistance developed within the beam element in response to the external loading (regardless of whether or not the response is material and geometric linear). Its evaluation involves numerical integration techniques, which will be discussed later.  Next is  introduced:  retmL  +  =  2 B 1 7fl  +  B 3 m  +  retmL  +  )dri 4 VB  re  =  IrLBdli  (2.41)  le  Then reiW =  where  reL  (2.42)  is the linear part of the internal force vector of the element and rede is the  non-linear part. The presence of the non-linear part is because of the geometric non linearity.  The work done by the inertia force in each element  lou TpiicIL? p°AüOu +-j p°NJONJ  pOJfr5()  +j  1O5  p0IO() +  pA17?OwIIx  (2.43) After substitution of the finite element version of the displacement field, for a typical element in the finite element:  IouTpudv = Ode Tmd  (2.44)  Ve  where  me is the element consistent mass matrix. For the higher order beam theory it takes  the form of:  27  me  =  [mi’]  [0]  [0]  [0]  [m2]  [m3]  [0]  [m2]  [m3]  (2.45)  where  1 m  =  16  C  23  me  1 Ndri I pAN  =  J  —  m  1 Ndn pIN =  dM• 32 Ndri =m efi p1 1 d  fl_i  33  me..  =  1 dM dM 1 + pAM, M p dq dq  (2.46)  jdii  The work done by the damping forces in each element can be similarly expressed in general nodal coordinates in the above way, i.e. IuTKdudv  (2.47)  8dCede  Ve  Then the work done by the external force in a typical element can be written as: fouTfdq = 6dre  (2.48)  le  wheref=jf,fI  and ret =  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 beam from the element generalized coordinates, connectivity transformations in the form of deCed  (2.50)  are employed to illustrate the procedure. In fact, the Direct Stiffness Method is often used in the practical application of the finite element method. Here d is the global vector of  28 generalized coordinates and C, is the connectivity matrix of bulean constants associated with each individual finite element under study.  Therefore, the virtual work statement, as expressed by Eq. (2.36), turns out to be as follows when expressed in the global generalized displacements  i3dT(CeTmeCed + CeTCeCe + CeTrt + CeTreCed)  =  ödT(CeTreet)  (2.51)  Since öd is arbitrary and kinematically admissible variations in the components of the generalized displacement vector, it can be cancelled out from the Eq. (2.51).  So the  governing differential equation of motion, as a result of the higher order beam theory and finite element discretization, is. md+cd+r =r  (2.52)  NE  where m  CeTmeCe is the global consistent mass matrix,  =  e=1 NE C  CeTCeCe is the global consistent damping matrix,  =  e=1  NE rIt  =  C’r + CTrJ1/jCd is the global internal force vector, and  e=1 NE  rt  CeTret is the global consistent mass matrix,  =  e=1  The governing equation (2.52) can be interpreted as external loads are equilibrated by a combination of inertial, damping, and internal forces. The material does not have to be elastic since the stress-strain relation is not specified during the deduction.  As observed, although the equations are spatially discretized, they are still a system of second order ordinary differential equations in time, hence continuous functions of time. Therefore, they are called a finite element semidiscretization. An explicit direct integration technique can be employed to discretize the system of equations in the time domain to  29 obtain a sequence of simultaneous algebraic equations [21]. This is discussed later in Section 2.5.  2.5 Modeling the Material Behaviour of the Beam  Before starting the solution procedure for the governing differential equations of motion presented 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 for each 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 Gauss  integration if the stresses are known. So, to determine r1t, the appropriate constitutive law for the material must first be derived so that the stresses can be calculated from the  strains. I  For ductile materials subjected to high intensity transient loading conditions, plasticity and strain-rate sensitivity will inevitably be a consideration. It is therefore desirable to incorporate elastic-plastic strain hardening and strain rate sensitive material behaviour into the material constitutive relation.  In the current study a bilinear elastic-plastic strain  hardening 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 the stress state is two dimensional, von Mises yield condition, the associated flow rule, and isotropic hardening are all employed to account for the plastic material behaviour of the beam. Ref. [22] gives a detailed discussion on the theory of plasticity, which is very helpful.  30 In plasticity, the following basic assumptions are based on experiments: At any stage of deformation, the strain tensor can be expressed as the sum of the  1.  plastic strain tensor  and elastic strain tensor —  where,  (p)  ,  (e)  253  4, by its mean g, is related to the stress tensor by Hooke’s law, (2.54)  =  2.  i.e.  The plastic volume change is negligible.  This is usually called the plastic  incompressibility, i.e. =0  (2.55)  To simpIi1r the material behaviour in the present problem, a bilinear elastic-plastic strain hardening model is used to describe the relation between stress and strain. Although the stress state in the present problem is no longer uniaxial, the bilinear model can still be applied to the effective stress and strain. The theory of plasticity is later used to determine the 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) with slope E yielding to a plastic strain hardening part (AB) with slope Et. That is, with a strain increment de, the stress increment is:  da = Eds  da = Ede or da = where ,as discussed before, is  for a<ay or da = Hde(1’)  =  +  for a>ay  (2.56)  and de is the total strain increment, c1e is  the elastic component of the strain increment and is recoverable, de’ is the plastic component and non-recoverable, and °y is the yield stress in uniaxial tests. The effect of  31  Fig. 2.3. An idealized bilinear model of the elastic-plastic material  32 strain rate on 0 y will be addressed later. Also  (2.57)  H= E  is the plastic modulus. If Et is set to zero, H will be zero, resulting in an elastic-perfectly plastic 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, and further plastic deformation is produced when the maximum previously applied stress is exceeded.  In addition to the bilinear elastic-plastic behaviour, the yield stress may now change in the current problem.  This is because the material is now subjected to high intensity  dynamic loading, which causes high strain rates in the material.  And it is shown  experimentally that, as the strain rate increases, the instantaneous yield stress of many materials increases from the quasi-static value. For some materials, this increase may be very significant when the strain rate is high. To account for this effect, the Cowper Symonds relation is used.  4i =  (2.58) +  where °dy and Gy are, respectively, the dynamic and static yield stresses. D and p are experimentally determined material strain rate parameters.  Summarizing, the modeling of material behaviour is based on the ideal bilinear elasticplastic linear strain hardening material models incorporated with strain rate sensitivity, as shown in the Fig. 2.4.  33  a  a 1 7 a  eo=o  e  Fig. 2.4. Idealization of elasticplastic linear strain hardening material behavior with strain-rate sensitivity From Ref [23J  34 To specify the constitutive relation, the yield condition must be known. Since the state of stress is multiaxial in the present problem, a yield fi.inction has to be used to identify different stages of material behaviour. Generally, the yield function depends on the state of stress and strain and on the history of loading, i.e.  (P)  41  ,K  where ic, the work hardening parameter, is a function of the plastic strain tensor  41.  Von  Mises yield condition is to be used in the current work. Introducing the equivalent or effective stress  as (2.60)  with s representing the deviatoric stress tensor.  In similar fashion, the equivalent or  effective plastic strain increment is  cI  =  j2%d$d4  The plastic work increment is dW  =  (2.61) Since the material under study is isotropic  with strain hardening, the yield function can be written as  f =f(o,H*())  (2.62)  where H* ((P)) is the plastic stress-strain relation between & and  Because von  Mises yield condition incorporated with isotropic strain hardening is employed, f may be written as  f =_H*())  (2.63)  After yield, the way plastic deformation further proceeds is governed by the associated flow rule, which states that the plastic strain increment must be proportional to the gradient of the yield function with respect to stress, i.e.  35  ôG,j  Ii  (2.64)  where ? is the proportional coefficient. The hardening rule governs the way the yield surface (the f=O surface in the stress space) changes in size and shape as the plastic deformation proceeds. For convenience, isotropic hardening (which assumes a uniform expansion of the initial yield surface) is used in the current analysis, although it has little experimental support. The hardening rule for the present material model is H*(1) =  (2.65)  b’  which is specified by the bilinear model.  As stated above, the yield thnction can be used to identifS’ different stages of material behaviour. With the yield condition well defined, the constitutive relation is ready to be defined.  For a state of stress corresponding tof(O, i.e. o< H, the material is elastic and there is no change in plastic deformation. The material stress-strain relationship is specified by Hooke’s law which is, in vectorial form for the current two dimensional stress state, as follows  dcr= [D]de  (2.66)  where [D] is the elastic constitutive matrix. This is applicable when the material is in the initial 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 plastic deformation occurs. This is the case when the material is loaded beyond the yield state or when the yielded material, after unloading, is reloaded beyond the previously applied stress state. Under these circumstances the stress-strain relationship is related to the  36 associated flow rule. Deduced from the theory of plasticity [23], the flow rule associated with von Mises yield condition in the form given in Eq. (2.63) is  .j 8 .C”) =(P)-i— ii aGU  (2.67)  In 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+A A=  1—T lôuJ  (2.70)  {Df1—T  Thus, 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.  The  incremental stress-strain relationship may be briefly put as  do= [D]d6  [D]tU ID I  1  =  °1  [D]—  [D]  H+A  (2.71)  where [D,, j is the elastic-plastic constitutive matrix for the total elastic-plastic strain increment.  The state of stress wheref0, i.e.  >  H, has no meaning and is thus inadmissible. To  summarize, Eq. (2.66) for yield functions less than 0 and Eq. (2.71) for yield functions equal to 0 comprise the elastic-plastic constitutive relation.  37  Since the calculation of the equation of motion is carried out step by step (which is discussed 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 state have to be scaled to the yield state of stress, and the pure elastic portion of the strain has to be subtracted, i.e.  —a 2 to=o 1 do- = =  where  (•)i  and  —  ( oi) +scal*du  ()2  scal=—= to ds = 82—81 (de)  =  (2.72)  (1— seal) *(d8)  refer, respectively, to the given quantity associated to the first or the  second state. After scaling, 01 is the stress state at yielding (which is a point on the yield surface) and de 1 is the corresponding total strain increment (which is composed of both elastic and plastic parts). Stresses can then be easily found using the plastic analysis for f=OinEq. (2.71).  Having established the constitutive relationship, the state of stress at any point of the beam may be determined given the strain. resultants i.e., 71,  , m, V 1 m . 1  Attention now turns to solving for stress  and V 2 along the beam axis, which leads to the evaluation of  the internal force vector. Because of the complicated distribution of these stress resultants in elastic-plastic analysis, the integral expressions in which they are specified are best integrated numerically. Among the numerical integration schemes, Gauss quadrature is employed in the current analysis because it is efficient, accurate, and suitable for the rectangular beam cross section. When written with respect to the natural coordinate, =  ,  the integration expressions for the stress resultants are as follows:  38  1 = -Jc*i m  (2.73)  =  where b and h are the width and height of the beam respectively.  Using Gaussian  quadrature [241: =  T7f() 1  (2.74)  ‘=  we have for the stress resultants:  2 m  2  =  = i=1  =  1=1  =V 1 +t  17t(j) =  (2.75)  1  which are enabled by the stress-strain relation presented before. Here, W 1 is the weight associated with the i-th depth-wise Gaussian evaluation point, and m=4. Four depth-wise Gauss points are sufficient to represent the non-linear stress distribution accurately as the material plastification extends through the beam cross section [9].  After the stress resultants are obtained from the stress by Gauss quadrature, the element internal force vector can be easily evaluated from the integration equation (2.40) and also  via Gauss integration. Written with respect to local coordinates 1+ ret” = J(nB using Gauss quadrature, gives:  +  3+ mB  +  flBde)dC  ,  =  —  1  (2.76)  39 reim  =  -  wj[n(C)B ( 1 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 as the finite element mesh is refined, three spanwise Gauss points are at least required within an element in order to accurately represent flBde and hence relm. Therefore n=3 [9]. With remt  known, r” can easily be assembled by using the Direct Stiffness method.  With all entries in the semidiscretized equation (2.52) known, the 2nd differential equation system in the time domain may be solved, and this is discussed in the next section.  2.6 Solution Procedure Over the Time Domain  To solve the dynamic differential equations over the time domain there are two types of methods: Modal Superposition and Direct Integration. For a blast loaded and non-linear problem, the Direct Integration method is favoured.  Among the Direct Integration  methods, explicit time integration techniques are well suited to treat the material non linearities. This method can be implemented easily and can handle very large problems with only modest computer storage requirements.  The only disadvantage is that the  required time step to ensure computation stability is small [21]. However, in the present blast loaded beam problem, the development of displacements and stresses are of interest  and a small At is necessary for accuracy. Hence the explicit time integration is used in the analysis. According to Ref. [9], the central difference is the most efficient method of time  40 integration when applied to problems involving impact or blast loading. Therefore, central difference is chosen to integrate the dynamic differential system of equations.  In the Finite Element Method, the spatial domain is discretized into finite elements. As a result, the finite element semidiscretization is:  [mjfi}+[c]{á}+{rmt} = frt} By using central difference, the temporal domain of the problem is discretized into equal time intervals of duration At. At each time interval a central difference approximation is used to replace the displacement time derivatives by differences of displacement {d):  f}  =  —  ) 1 fd}_  fJ}  =  1 ._--({d}÷  —  2fd} + {d}_ ) 1  (2.78)  With the central difference method, the recurrence relationship [211 is  m+  —-  c][d}÷  =  {rt  }  —  frmt} + ._!j..[m](2fd}  —  )+ 1 {d}_ (2.79)  whereby the unknown displacements at t+At are given in terms of known quantities at  time t. Since this is a two step method, a special starting procedure is required. The starting conditions are as follows:  =  =  As  0 fd}  —  0 +.-fi} &{a} 0  [m]_1({relt}o  observed from the  recurrence  —  {rh1t}  —  [c]fii} ) 0  formulations,  central  (2.80) difference  clearly is  computationally simple, since the unknown displacement vector occurring at any time step  41  is directly obtainable by performing simple arithmetic operations on matrices of known value.  If [m] and [c] are diagonal, then the equation system is uncoupled and the displacement vector can be obtained without solving simultaneous equations, thus making the computation even simpler. Actually, according to Ref. [21], explicit integration is usually more accurate with lumped mass matrices than consistent mass matrices. Therefore, with central 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: y’ 1 PA 0  io  /  pAle 3 /24  rLmeJ_ i_  pAle/ /2 34  0  105 PAl/ pAle / 3 /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 and accuracy of the computation.  It is to be noted that central difference is conditionally stable and requires At such that:  42 (2.83)  AtAtcr  even when the structure is restricted to being linear. Here, Atcr is the critical time step to ensure computational stability according to Courant’s criterion =  2  (2.84)  Atcr  max 0 If At is exceeded in time integration, the computations will be unstable resulting in erroneous time-histories as they grow unbounded.  is the highest natural frequency  occurring 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 nonlinearity does the opposite, an estimate of the critical time step is usually obtained by reducing the corresponding linear result by 10-20% [9]. The present analysis can be further simplified because omax for the assembled finite element model is bounded by the maximum natural frequency of the smallest constituent unassembled and unsupported element, and this frequency can usually be computed easily [21]. For an Euler-Bernouffi beam, the critical time steps are as follows: for consistent mass matrix: Atcr  =  &cr  le/(q5cL)  =  if le  l/(10JJcLrG)  1oq’Yr otherwise  for diagonal mass matrix Atcrle/CL Itcr =  where  CL  if le4. rG 5  l/(4...JcLrG)  otherwise  (2.85)  is the speed at which longitudinal waves propagate along the beam axis  =  VP  rT.is the radius of gyration of the beam cross section. .  and rG  .  In the higher order beam  =  theory, shear deformation is taken into account and the beam is not as stiff as an Euler Bernouffi beam, which means the highest natural frequency occurring in the higher order  43 beam will be lower. Therefore using the stability bound for Euler-Bernoulli beam element is conservative, which is desirable.  As observed from the above stability analysis, lumped mass matrix clearly provides for larger stable time steps in addition to providing uncoupled equations and generally more accurate results than the consistent mass matrix in the central difference method. It is also to be noted that the critical time step in time integration is governed by the size of the smallest element in the mesh, and the refining finite element mesh will require a more restrictive critical time step.  Therefore, extreme mesh refinement in any part of the  domain is undesirable.  According to Ref. [21], the stability requirement is governing for the explicit time integration method. For systems of finite element equations, excellent accuracy can be arrived at by using a time step just under the stability limit. So in the present central difference method, a time step that satisfies stability criteria can usually guarantee accuracy satisfactorily.  With the discussion about the integration scheme completed, the finite element semidiscretization may be turned into a system of uncoupled algebraic equations to solve for the displacement vector. Moreover, by using the central difference, the solution in the time domain may be used to obtain the dynamic transient response of the structure.  2.7 Summary  This chapter presented the formulation of the numerical scheme to simulate the blast loaded beam problem.  Emphasis has been laid on the kinematics and the governing  44 equation of motion of the higher order beam theory, and elastic-plastic modeling of the material with strain rate sensitivity.  With the whole formulation complete, the elastic-plastic transient response (including the non uniform shear effect) of the blast loaded beam may be numerically simulated. This provides the basic data for the failure analysis later.  45  Chapter 3 Failure Analysis of the Problem  3.1 Introduction  Using the formulation discussed in the previous chapter, the transient response of a slender ductile beam with geometric and material non-linearities under blast loading conditions may be numerically simulated. With the displacement, stress and strain fields of the beam available, only the failure criteria is needed to predict failure when the beam is subjected 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 previous research work by Fagnan [10], several versions of interaction failure criterion are developed. Then, the failure condition using the theory of plastic work density is also formulated. 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. discussion.  The results are presented in the next chapter for  46  3.2 Interaction Failure Criterion  Menkes and Opat [6] conducted a series of experiments on the dynamic plastic response and failure of fi.illy damped metal beams which were subjected to explosive loading conditions. At low impulsive loading conditions, large permanent ductile deformations of the beam were observed.  At intermediate loading intensities tearing rupture was  encountered, which was caused mainly by excessive tensile strains developed at the supports of an axially restrained beam. At still higher loading intensities, a more localized failure occurs owing primarily to the influence of large transverse shear forces. Between the threshold values for the two failure modes, the beams failed with a mixed tensiletearing and transverse-shear mode. This indicates that it is possible that both the tensiletearing part and transverse-shear part contribute to the failure of blast loaded ductile beams.  This is the idea behind the interaction failure criteria.  In this criterion, the  contribution of the tensile-tearing part is specified by a tensile strain ratio and the extent of the influence of the transverse shearing part is taken care of by a shear stress ratio. Both parts are then included in the failure conditions to predict the failure of the beam.  3.2.1 Calculation of Influence of Tensile Tearing  In the interaction failure model, a tensile strain ratio is used to describe the contribution of the tensile-tearing behaviour of the beam to the failure of beam.  Since the tearing  action of an axially restrained beam is most severe at the support, the maximum strain of the beam at the support is needed in the tensile strain ratio expression.  47 The elastic-plastic finite element analysis F.ENTAB version 2 [13,25], which is programmed according to the procedure presented in Chapter 2, is able to give the displacement, strain, and stress field of the present beam problem. A typical deformed profile, obtained by FENTAB, of a clamped beam subjected to a uniform lateral load is shown in Fig. 3.1. However, as seen from the figure, the program is not able to accurately model a plastic hinge at the support of the beam. A modified approximate shape of the beam which accounts for the plastic hinge is therefore needed to predict the maximum tensile strain at the support.  As mentioned in Ref. [26], Onat and Shield suggested that there exists a triangular plastic region at the support of the beam under impact conditions. This suggestion is encouraged by the appearance of deformed aluminum specimens in experiments [27]. Thus, while using rigid plastic theory, Jones [7,8] modelled the clamped beam under impulsive loading as two rigid parts joined by a plastic hinge at mid-span and two plastic hinges at the supports.  Based on this approach, Fagnan [10,27] derived a model  consisting of plastic hinges only at the supports and the rest of the beam as modelled by the Finite Element Analysis. This approximate shape of the beam is adopted in the present work 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 the support, it is assumed that the elastic response is negligible compared to the plastic deformation.  With the modified deflected shape of the beam established, the maximum total strain at the supports can be developed from the beam shape model. In the present problem of a fully clamped beam, the maximum total strain 6 tot at the supports can be expressed as follows tot 6  bend +8ial 8  (3.1)  48  —  L  —  I —I 9  Fig.3. 1 Typical deflection profile as determined by FENTAB. From Ref.[1O]  49  where  ej  is the axial strain at the neutral axis at the support and 8 nd is the bending  strain at the support.  Here, for a slender beam, the bending strain may be expressed as  bend 8  =  Kh  (3.2)  where h is the depth of the beam and K is the curvature of the plastic hinge. For a plastic hinge, the curvature K 15 defined as (3.3) ‘ph  where 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 a plastic hinge, approximation is made for the plastic hinge length 1 ph A slip line field analysis [26] can be used to obtain the shape of the plastic hinge across the depth of the beam. It reveals that a plastic hinge in a beam has the shape shown in Fig. 3.2(a) when the beam experiences pure bending behaviour and has infinitesimal transverse displacements. As the mid-span displacement increases, the membrane force increases and the axial stretching of the plastic hinge occurs until it takes the form in Fig. 3.2(b) with the onset of a membrane response (when the mid-span transverse displacement equals the beam thickness h for rigid plastic analysis on a fully clamped beam). As observed from the figure, it is evident that the length of the plastic zone on the upper surface of a fhlly clamped 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 an approximation for the hinge length “ph is lph”th  (3.4)  50  l=h  (Q)  F  lt2h -  (b)  Fig. 3.2 Plastic hinge in a fblly clamped beam: (a) pure bending state; membrane state. From Ref.[12]  (b) onset of  51 where a is the plastic hinge parameter whose value is chosen by the user. According to the study done by Nonaka on the hinge length [26], for rigid plastic beams of rectangular cross section, 1 a  2 for maximum transverse displacement between 0 and h. After that  the beam is regarded as being in the membrane state, the plastic deformation is assumed to develop all over the beam, and the plastic hinge length is assumed to be: ‘ph  =  where  L is the half length of the clamped beam. Jones [11,121 suggested that the plastic hinge parameter change with the impulse according to an empirical formula. In the present elastic-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 the finite element procedure is not able to model a plastic hinge at the support.  The  assumption is made that the rotation of the plastic hinge is equal to the maximum rotation occurring in the beam, which is usually a short distance from the clamped support as seen from Fig. 3.1. To determine the maximum rotation of the beam, first the generalized nodal coordinates, as given by Eq. (2.30)  de  I  1 &w  r 2 w 8 1  at all the finite element nodes can be examined, and the maximum nodal rotation value  ox  can be easily selected. Second, to be more accurate, interpolation of the nodal value can be made along the length of each element using the shape functions to predict the maximum 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 as  4 , 2 wfO,0,Mi,M , 3 } 0,0,M M de  (3.5)  1 M 2 , 3 M are cubic Hermitian polynomials as specified in Section 2.4. The where 4 rotation of the beam can then be expressed by the first derivative of w with respect to x,  52  that is with respect to 11 (local element coordinate) within each element  —10 Setting  —i-  =  0  2 1 8M 8M  0 0  3 8M 8M 4  ‘d  36  0 and back substituting, the location and magnitude of the maximum  rotation within each element can be determined. When written with respect to the local 2’ —1, they are as follows. normalized coordinate  c  (ãw 1  ‘2  —  Cmax  —  6(w1_w2)+3le(-1+-)  = —  1)W1 +(3C2  —2C— i)  2 +. L+(i_c2) (3C2 4  +2C—  ()max  (3.7) The interpolation procedure can be carried out for each element so that the maximum rotation  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:  Ebend  0 2ci  (3.8)  Attention now turns to the evaluation of the axial strain along the neutral axis at the support of the beam.  This axial stretch can be obtained by directly integrating the  deformed profile as determined by finite element analysis of the beam element closest to the 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:  53  I  10  1=  ldx)  0  After Taylor expression gives 1 l0(f12  +•J  llo  (3.10)  from which the axial strain is obtained 1 axial = 8  (3.11)  2l  This is the formula that is used to calculate axial strain at the boundary. Noticing the finite element approximation of the displacement field (3.5) within each element, for the first element: 1 w=fO,0,M , 2 , 3 } 4 o,0,M M d  (3.12)  where d 1 is the displacement vector associated with the first element closest to the boundary. Then  ax  in the first element is the same as  100  1 8M 8M 2  00  3 ôM  (3.13)  Next T  2 8M  110  1 4 oM Eial_JdlT0,0 ‘o’aq’’’aq’aqJ 2l (  2 8M  0  (3.14) recognizing B  —  {o 0  1 8M 8M 2  0  3 aM aM 4  1T1 0  1 8M ôM 2  0 0  3 ?jM 8M 4  1  —  Expressed with respect to the local nondimensionalized coordinate axial strain at the support is ax,al _!‘J•diBdidc 6 —  (3.15)  ,  C  =  —1, the IC  54  Choosing three spanwise Gauss points, which is required by the exact evaluation of B, the integration by Gauss Quadrature is evaluated. In this way, the axial strain at the boundary  is determined.  The total maximum strain at the support is now available by summing the bending and axial 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 as the contribution of the tensile tearing behaviour to the damage of the beam.  3.2.2 Measure of the Influence of Transverse Shear  Measuring the influence of the transverse shear is accomplished by using a shear stress ratio. According to the experiments, this influence is most significant at the support, and  metal beams are sheared off at the support during mode ifi failure.  Therefore, to  completely consider the influence of transverse shear, the shear stress at the support has to be evaluated. In the present finite element analysis the higher order beam theory (which considers non uniform shear distribution across the beam depth) is employed, allowing analysis of the shear force and shear stress of the beam. As observed from the shear force distribution of the beam, however, the shear stress and resultants are more accurate in the middle of the element than elsewhere. Shown in Fig. 3.3, the shear stress distribution fluctuates in the wave length of about one element size. This is probably because of the  approximate nature of the finite element method. To overcome this detriment of finite  55  2  .0  o  •1—  i Z  0  0  1  2  3  Distance from the support (in.)  4  5  4  5  1000  800  600 .0 C)  0  —200 0  1  2  3  Distance from the support (in.)  Fig. 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.  56  element analysis, two options are available. One is to obtain the shear stress by evaluating the reaction force at the support. The other is to use the shear stress at the middle of the boundary element to approximate the shear stress at the support. The first method shall be 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 is required. 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 the external 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 the equilibrium can give the end reaction in this situation. The above discussion can be put into the following expression. L  L  L  R + jq(x)dx—Jm(x)dx—Jc(x)i4’dx= 0 0  0  (3.17)  0  where R is the wall reaction in the vertical direction, q(x) is the vertical component of the external distributed load, m(x) is the sectional mass density of the beam with a unit of mass 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 neutral axis of the beam.  As a result of applying the finite element method, the continuous beam is discretized into an equivalent multi-degree of freedom system with the degrees of freedom {d}, associated mass [ml, damping [c], and applied load {rt). Since lumped mass and damping are used, the equivalent multi-degree system reduces to a lumped mass system with mass and damping associated to its own degree of freedom at each mode. equilibrium equation turns into:  Therefore, the  57  q.  R +  nodes  131 m  —  nodes  13 =0 c 1 *  —  (3.18)  nodes  where qj is the effective external force on each node, mj 3 and c 13 are the mass and damping associated with vertical degree of freedom at each node, respectively, and  and  1 are the vertical acceleration and velocity at each node, respectively, which can be w obtained from the transient finite element analysis. In this way, R can be found at each time step.  With the shear reaction at the wall known, the shear stress is easily derived. For a rectangular beam at an elastic stage, the maximum shear stress is at the centroid of the = 1.5Rw/ beam section with a value of where A is the area of the cross section. ,  However, as the central fibers of the beam yields, the shear stress distribution across the depth of the beam will become more uniform. In addition, it is assumed that the shear dominant failure will not occur until all the fibers in the beam cross section reach the ultimate shear strength. Therefore, the shear stress in the present analysis is treated as being distributed uniformly across the beam depth at rupture of the beam, and can be obtained from the wall reaction as (3.19)  Finally, to measure the influence of the transverse shear force, the normalized ratio  (3.20) tult  is employed, where tult is the ultimate shear strength of the beam when the beam is in the pure shear state. As with the tensile strain ratio, the transverse shear ratio reflects the amount of shear damage.  58 3.2.3 Interaction Failure Criterion  The previous two sections discussed how to evaluate the influence of tensile tearing and  transverse shear on the failure of the beam. As stated ahead, a tensile strain ratio and a shear stress ratio are employed to measure the influences.  These two ratios are also  reflections of the damage extent caused by, respectively, tearing and shearing.  With the measure of the two factors causing failure defined, the failure conditions may be specified. The failure of the structure is related to tearing and shearing actions of the beam. 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 be  ‘‘  jEj  t)  comprehended to be a general equivalent degree of the structure damage.  As it is  expected failure will occur when the extent of damage is 100%, i.e. 1, it is then reasonable to assume that failure of the beam occurs when  (3.22)  tft)  Ef  For the failure thnctionf two functions are assumed:  tmax  1)  I  Ej  t,jt  )  +  =  tuft  Sf  ‘‘  tmax  (3.23)  2,  2) Sf  taft)  8.1)  taft  59  which are the linear interaction failure criterion (LIC) and the quadratic interaction failure criterion (QIC).  3.2.4 Another Two Interaction Failure Criterion  As mentioned is Section 3.2.2, the shear stress is obtainable from the finite element  analysis because the shear effect is included in the analysis procedure. However, the shear stress is not accurate at the support because of the approximate nature of the finite element method. Since sufficiently accurate shear stresses at the mid point of the element are obtainable, it is possible to approximate the shear stress at the support using the shear stress at the middle of the boundary element. In this way, as the finite element mesh near the 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 center fibers in the section yield. It is also assumed the shear dominant failure does not occur until all the fibers in the section reach the ultimate shear strength. In addition, there might be possible fluctuations in shear distribution across the beam depth in the results of finite element analysis. Thus the average stress is used JtcL4 tavg  where V  =  (3.24)  ftdA, the first order shear resultant force. A  Substituting tavg for  in the two failure conditions (3.23) in the previous section,  another two interaction failure criteria are derived  60 1)  f1L, 1  Cf  tayg  +  =  tuft)  tjjf  f 6  (3.25)  2,..  I  2  +‘‘  2)  taft)  j Cf  I 8 f)  taft  3.3 Failure Criterion Using Plastic Work Density Theory  The work density caused by the stresses at each point of a structure consists of two parts: 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 of plastic 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 stress components, and cr is the critical value of plastic work density. It can be taken as  4?) =  where  ad(C,é)  f  (3.27)  Od(C,6)dC’  is the dynamic engineering stress, which is obtained from a dynamic  uniaxial tensile test for a given é, and where s is the plastic strain at rupture. Since elastic strain is far less than plastic strain at failure, plastic strain at failure is taken to be equal to the rupture strain, i.e.  4?)  =  (3.28)  where e is the engineering rupture strain. Actually ad and  magnitude of the strain rate, in general. However,  Cf is  both could depend on the  assumed to be independent of ê in  the present analysis. The Cowper-Symonds relation is used to take strain rate sensitivity  61  into account in a (E,). Since a bilinear material is assumed for the elastic-plastic strain hardening material, c 0 r can be expressed as follows:  cr dyef  (3.29)  . 1 +--Es  where Et is the tangent modulus of the material and aê can be obtained from the Cowper Symonds relation discussed in Chapter 2  a,=a ‘  (flY  1+1— D  As to the density of plastic work, S can be similarly expressed as E(p) =  fajcic””  (3.30)  where ad and e(P) are true dynamic stress and plastic strain, respectively, in an uniaxial case. In an actual structure, they are assumed to be equal to the equivalent stress and strain. In the present finite element analysis the stress and strain state at Gauss points are available at each time step. Thus, 0 d and zSE can be obtained at each time step. S, the density of plastic work at a point, can therefore be evaluated by trapezoidal integration. (p) 8 =  (3.31)  where aj 1 and a 1 are the stresses at i-1-th and i-th time step, respectively, and zSz is the plastic strain increment at i-th time step.  When the density of plastic work is greater than the critical value, a crack is initiated at that point. However, after a crack is initiated but before severance occurs, the structural  62 member might continue to support loads. The beam is therefore assumed not to fail until the whole section fails, that is, until the plastic work density in the whole section reaches the critical value.  cl=.clcr  (3.32)  where 2 is the sectional plastic work density, whose unit is work per length. cr is the critical value.  Because the plastic work density at each Gauss point is known, the sectional plastic work density can be calculated. The sectional plastic work density is the plastic work density in the whole section, i.e.  c=JJ&L4  (3.33)  A  For a rectangular beam with width b and height h,  f1=b8dz  (3.34)  Expressed in terms of local nondimensionalized coordinates, bh 1  where  =  .  (3.35)  Then the sectional plastic work density can be evaluated by uss  quadrature. Four Gauss points are used in the Gauss integration. Therefore, bh  ‘  2 i1  (3.36)  where 9. is the plastic work density at i-th Gauss point. W() is the weight factor associated with i-th Gauss point.  63  Similarly, the critical value for the sectional plastic work density can be written as  cr  =JJcr’4  (3.37)  A  When expressed in terms of local nondimensional coordinates,  bh 1 cr=fcr’1 2 1  (3.38)  Because at each point of the section the strain rate might be different, the dynamic yield stress at each point might be different and so is the critical plastic work density. Therefore to evaluate ncr’ Gauss integration is employed. Applying Gauss quadrature,  2cr  where  bh 2  ‘  (3.39)  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 will fail when, at one section, (3.40)  This comparison between 2 and cr can be carried out at each spanwise Gauss point along the beam to decide whether the beam fails or not.  64  3.4 Post Failure Analysis  With the finite element procedure established in Chapter 2, and the failure conditions discussed in the present chapter, the transient response and failure of the blast loaded beam can be simulated. However, it should be noted that the transverse displacement of the beam at the instant of rupture is usually different from the permanent transverse displacement obtained from the experiments. This is because the beam has a residual energy at severance when all of the initial energy input has not been absorbed in plastic deformations completely before failure, which is usually the case. The residual energy at rupture is often significant. With so much energy at failure, the beam is bound to deform until a stable state of the beam is reached.  As the beam continues to deform after  severance, part of the residual energy is absorbed in fhrther plastic flow until it reaches the steady state where no more plastic deformation occurs. At the steady state, the beam moves with associated rigid body kinetic energy and vibrates with relatively small amounts of energy exchange between elastic energy and kinetic energy, simultaneously. The midspan transverse deformation at the steady state, which is the average of the maximum and minimum 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-free beam with initial motion. After severance of the beam at failure, the beam is blown off and becomes a free flying beam. To numerically simulate this, the boundary which was clamped before is now set free. Further temporal progression of finite element analysis for a free-free beam can give information on the motion.  After incorporating the post failure analysis into the program, it is now possible to simulate 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 experimental  65  data obtained by Menkes and Opat[6] and the discussion of the numerical model are presented in Chapter 4.  3.5 Summary  Failure conditions are developed in this chapter. The interaction failure criteria and the failure criterion using plastic work density theory both consider all contributions from tensile and shear action of the beam.  In addition, to fully simulate the response of a blast loaded beam, and to check with the experimental results properly, post failure analysis is discussed.  This concludes the theoretical modelling of the problem.  With the theoretical base  presented in Chapters 2 and 3, the whole response process of the beam may now be simulated numerically.  66  Chapter 4 Numerical Example & Discussions  4.1 Introduction  With the finite element procedure discussed in Chapter 2, and failure conditions and post failure analysis established in Chapter 3, it is possible to obtain the transient response of the blast loaded beam, to predict the failure, and to obtain the mid-span transverse deformation in the steady state. However, the efficiency and accuracy of the prediction is still unknown.  Menkes and Opat [6] conducted a series of experiments on 606 1-T6 clamped aluminum beams subjected to a surface explosive charge. These experiments are used as the basis for comparison.  To check the efficiency and accuracy of the numerical model, finite  element numerical simulations of these experiments are carried out. Calculation results are then compared to the experimental data, giving encouraging correlations. Furthermore, some insights into the response of blast loaded beams are obtained through the finite element analysis.  67  4.2 Experiments  Menices and Opat [6] conducted a series of experiments on clamped aluminum beams subjected to blast loads in which the uniformly distributed impulsive load is provided by sheet explosives. The experimental configuration is shown in Fig. 4.1. As shown in the configuration, the sheet explosive is cemented to a neoprene buffer which is bonded to the top 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, and uniformly distributed impulse loads can be assumed. The intensity of the impulse was varied by using different combinations of sheet explosive of different thickness. The beam specimens 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 different damage modes can be identified for the beams with rectangular cross sections, as shown in Fig 1.1: I Large inelastic deformation of the entire beam, II Tearing (tensile failure) of the beam material at the extreme fibers at the supports, and ifi 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 acquired excessive permanent transverse displacements without any material failure. The damage  severity can be measured by the residual central deflection. At mode II the beams failed due to dominantly tearing of the beam material at the supports. The threshold is taken as the 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 severed central section decreases (which can only be measured from photographs) until mode ifi is  68  —  EXPLOSIVE —  NEOPRENE BEAM  iI  2L  Fig. 4.1. Schematic drawing of experimental configuration. From [6,10]  69 reached. This is shown in Fig. 4.2. In mode ifi, failure of the beam occurs because the influence of the transverse shear forces dominates the response and causes failure of the material at the support. The threshold is taken as the impulse intensity which first causes no significant deformation in the severed central section. Plastic deformation of the beam becomes localized near the supports and no appreciable plastic deformation is present over most of the beam at this stage. The permanent deformation of the severed central section remains almost constant. When the beam is subjected to loads which lie between the mode II and mode ifi thresholds, failures which involve both the tearing and shearing modes are observed.  4.3 Modelling of the Problem  The modified finite element analysis FENTAB is used to numerically simulate the blast loaded beam experiments conducted by Menkes and Opat  .  For simplicity, the actual  explosive loading condition can be modelled as a rectangular pressure pulse, with a constant 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. The validity 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 transverse displacement is insensitive to this dynamic pressure ratio. This suggests that the details of the pressure pulse are not important [1,81. According to Fagnan [10], the duration of the rectangular pressure pulse is chosen to be 60 sec.  Since the load and the clamped beam are symmetric, only a half span of the beam is studied. To represent the half span in the calculation, several finite element meshes have been used. Ten elements of equal lengths are used most often. When a finer grid is  70  KIAPS 10.9 17.8 25.6 28.7 35.1 39.6  42.9 46.0 50.7 52.9 6 .6  Fig.4.2. Experiment results for series of 6061-T6 beams (0.25x1x8-in. beam) From [6]  71  needed, four different meshes were used: (i) sixteen elements of same size, (ii) twenty elements of equal length, (iii) twelve elements with the first two elements of the ten equal element mesh subdivided into four equal length elements, and (iv) fourteen elements with the first two elements of the twelve element mesh fi.irther subdivided into four equal length elements.  The material of the beam is modelled as elastic perfectly plastic, i.e. Et=0. The reason for the modelling is that (according to Jones [11]) during the dynamic response of the beam the strain rate will decrease with time, resulting in a stress-strain curve which approaches perfect plasticity. Other properties of the 6061 T6 aluminum as quoted by  Fagnan [101 are p=2.56x10 4 /in 41 2 lbs , 6 psi. Since 6061 0o 6 , 00 psi, and E=10.5x10 T6 aluminum is essentially strain rate insensitive at the strain rates encountered in the experiments, i.e.  =  o,. In the current study, damping is set to zero because the most  severe condition of deformation is desirable for the analysis.  To predict the failure of the beam, the ultimate data of the material is needed. The maximum allowable strain is chosen to be 0.17, which is the percentage elongation at fracture of a static uniaxial tensile specimen made of 6061 T6 aluminum. Menkes and Opat [6] observed from the experiments that the dynamic shear strength appears to be much higher than the static one.  After considering the von Mises yield criterion and  ultimate tensile strength of the material [10], the ultimate shear strength is chosen to be 30,000 psi.  With the experimental conditions modelled, results can be obtained from the numerical simulation, which are discussed in the next section.  72  4.4 Numerical Results and Discussion  4.4.1 Numerical Simulation of the Response Process  After modelling the problem of a blast loaded beam as described in Section 4.3, using the finite element procedure which is presented in Chapter 2, incorporated with failure conditions discussed in Chapter 3, it is now possible to numerically simulate the whole response process of the beam subjected to blast loading conditions. An example of this is discussed below.  Fig. 4.3 shows a series of proffles of one half of a blast loaded beam. The beam is 0.375 in. high, 1 in. wide, and 8 in. long. The impulse is 52.9 kLaps (1 ktap 1 ktap  =  =  0.014 lbf-s/in 2 or  100 2 N-s/rn ) . As observed from Fig. 4.3, the part of the beam close to the  boundary undergoes severe deformation while the middle part of the beam experiences little deflection in the early response.  As time progresses, the deformation gradually  extends to the middle of the beam. When certain failure conditions are satisfied, the beam breaks away from the supports, becoming a free beam.  While flying away from the  supports, the beam experiences further plastic deformation, resulting in more energy being absorbed 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 is the difference of the displacement of the mid-span node and the boundary node after rupture. 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 the highest peak, more energy is absorbed by plastic work, and then the beam reaches the  73 3:5  3 C  2.5  C 0  E  0  0 0  1.5  I  0.5  0 0  3  Distace from origin (in.)  4  0.05  —  0.04  C C  I0 0.02 0 C 0  0.01  0 0  Dista,ce from origin (in.)  025  —  02  C  = 0.15 0.  0  :  0.1  0 C  0  0.05  0 0  1  2  ..  Distance from origin (in.)  3  4  Fig. 4.3. Response profiles of one half of a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 52.9 ktaps.  74 0.7  0.6 C — 0.5  C 0  E 0  0.4  0  0.3 0  I: 0  0  1  DistLice from origin Qn.)  45  4 C C 0  E  0  3.5  —  0. 0 0 0  3  0 C 0  I 2.5  -  550 microseconds, close to 1St highest peak of central deformation  u=0 40659 in. 2  1  0  5  Distnce from origin fin.)  8  7.5  —  C C 0  E  0 0  7—  a  0. 0 0 0 0 >  6.5  -  0,  C  a  I-  8 -  1050 microseconds, close to 1st lowest peak of central deformation  u=0.24904 in.  5.5  I  0  1  .2  ..  Distance from ongin fin.)  3  4  Fig. 4.3. Response profiles of one halt of a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 52.9 ktaps. (continued)  5  75  2.5  2  C C  o  1.5  .  1  (3 0.5  0 0  1  2  3  4  5  6  Time (millisec.)  7  8  9  10  Fig. 4.4. Time history of central deformation for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 52.9 ktaps. The failure is predicted by the quadratic interaction failure criterion (QIC).  11  76 steady state. The permanent deflection, which is taken to be the average of the maximum  and minimum central deformation Aw in the steady state, is 1.474 in. Obviously, after rupture, the beam experiences further severe deformation because of residual energy at the rupture of the beam, which cannot be neglected. In this case, because of rupture of the  beam, the permanent deflection is less than the deflection when the beam is considered to have indefinite tensile strain.  Further study on the steady state of motion reveals that the period of the free beam broken away from the supports, as obtained from numerical results, is 820 .tsec. For a uniform Bernouffi-Euler beam with both ends free, the circular frequency of first non rigid body mode is [28]: 2 CD  =  (4.73 1)2.Vi-  (4.1)  which turns out to be 820 .tsec. Since the mid-span deflection time history of the beam without rupture obtained from the higher order analysis is almost the same as from the Bernoulli-Euler beam analysis, as shown in Fig. 4.5, there should be no significant difference between the natural periods of the two beams. The natural frequency of a  higher order beam with both ends free should be around 820 p.sec, which agrees with the result from the numerical simulation.  Compared to analysis without considering the shear effect, the time history of mid-span deformation 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 early stage of response. Shown in Fig. 4.7 are some proffles during the early stage of response of one half of a 0.3 75 in. x 1 in. x 8 in. beam subjected to an impulse of 42.9 ktaps. Since the 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 more  77  0.8  0.7  0.6 C  E 0 C 0  0.2  0.1  0 0  0.1  0.2  Results from FENTAB  0.4  0.3  Time (millisec.) —--  0.5  0.6  0.7  Results from FENTAB v.2  Fig. 4.5. Comparison of time histories of midspan deformation obtained 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 17.8 ktaps.  78  0.5  0.4  0.2 0 Co  0.1  0 0  0.01  0.02  0.03  0.04  0.05  Time (millisec.) __  Results from FENTAB  _—  0.06  0.07  0.08  0.09  Results from FENTAB v.2  Fig. 4.6. Comparison of time histories of midspan deformation during early stage of response obtained 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.  79 1.5  0.5  0  .....  Results from FENTAB  _..,..  Results from FENTAB v.2  0.04  0.03  C C  :0.02 C) 0  = 0  0 0.01  0  _..  Results from FENTAB  ...  Results from FENTAB v.2  0.25  0.2  0.15 C  0  0 0  0.1 0  0.05  0  Dista i 3 oe from origin (in.) ...  Results from FENTAB  ..  Results from FENTAB v.2  Fig. 4.7. Comparison of beam profiles during the early stage of response obtained from FENTAB and FENTAB v.2.  80 severely in the early response.  After the different responses in the early time, the  difference of the profiles close to the boundary between the two analyses decreases as time progresses. Finally, if the beam can sustain the impulsive load, the difference will be very small.  4.4.2 Discussion of Failure Conditions  From the above discussion, the response processes of a beam subjected to different loading conditions are now available. However, different failure conditions give different predictions of time of failure, or even different predictions of whether the beam can sustain the impulse or not.  That is, different failure criteria predicts different processes of  response. Therefore the accuracy of different failure conditions in predicting the failure  needs to be addressed. This is discussed in this section.  1. Interaction Failure Criteria  After analyzing the 0.375 in. x 1 in. x 8 in. beam subjected to different impulsive loads by using the interaction failure criteria, the mid-span permanent deformations (which are the differences between the displacements of the mid-span and the supports) are plotted against the impulses in Fig. 4.8. Experimental results of mid-span permanent deformation is available from Ref. [6] for the beams which did not fail. Results for the broken beams are digitized from the figure in Shen and Jones [11], which they obtained from measuring the photographs in Menkes and Opat’s experimental report. For comparison, these results are 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 permanent deformation increases as the impulse increases. Then, after a certain impulse value, the  81  2.5  C 0 (U  1  0 0  10  ±  20  40 30 Impulse (ktaps)  50  QIC (10 elements)  LIC (10 elements)  QIC (20 elements)  QIC (16 elements)  Mode I Exp. Results  60  70  x Shen & Jones’ (Mode II Results)  Fig. 4.8. Midspan permanent deformations vs. impulses obtained from a 0.375 in. by I in. by 8  in. beam using the quadratic interaction failure criterion (QIC) and the linear interaction failure criterion (LIC).  82 mid-span permanent deformation starts to decrease until the beam reaches mode ifi  failure, where the deformation is almost constant. However, this criterion predicts mode II tensile rupture of the beam far too early at 21.8 ktaps, as opposed to 46 ktaps reported by Menkes and Opat. It also underestimates the mid-span permanent deformation of the broken beams.  Therefore, the linear interaction failure criterion cannot predict the  response of the blast loaded beam as accurately.  Using the quadratic interaction failure condition, results of failure analyses employing three different finite element meshes (10, 16 and 20 equal elements) are presented in the figure. The curves also predict the trend of mid-span deformation correctly. It fits the experimental data very well and is on the conservative side. A detailed look at the failure of the beams reveals that the first rupture of the beam occurs at 289 IJ.sec. According to Menkes 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, the quadratic interaction failure criteria still gives a fairly early mode II rupture. Experimental data of Menkes and Opat shows that the mode II tensile rupture occurs at 46 ktaps for the 0.3 75 in. x 1 in. x 8 in. beam, while the quadratic interaction failure criteria predicts the mode II threshold to be 28.7 ktaps.  To predict the mode ifi failure, the tensile strain ratios ratios  and the shear stress  are plotted versus impulses. The results of three quadratic interaction failure  analyses employing different grids are presented in Fig. 4.9. As seen from the figure, the tensile strain ratios are greater than the shear stress ratios and they are almost constant during the mode II failure. Then, at a certain impulse, significant increases or decreases appear in the curves of tensile strain ratio or the shear stress ratios. The ratios of shear  83  1 tensile strain ratio  0.9  0.8 U,  0  0.7  -  0.6 shear stress ratio  0.5 I  0.4  20  I  30 _._  __-  tensile ratio (10) shear ratio (10)  40 Impulse (ktaps) _-•-_  __  tensile ratio (16) shear ratio (16)  I  I  50  60 —A-__  70  tensile ratio (20) shear ratio (20)  Fig. 4.9. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion (QIC).  84 stress exceed the ratios of tensile strain and start to be the controlling factor in the failure of the blast loaded beam. The impulse immediately after the big increase or decrease is specified to be the threshold of mode III. For the 0.375 in x 1 in x 8 in. beam, it is 50.7 ktaps,  while the threshold for mode ifi failure is reported to be 64 ktaps.  This  discrepancy is probably because of different definitions of mode ifi threshold. The mode ifi threshold in the present analysis means a starting point where the shear stress ratio  starts to dominate, while the mode ifi threshold in Menkes and Opat’s experiment is taken as the impulse intensity which first causes pure shear failure with no significant deformation of the severed central section. The strain and stress ratio chart for analysis using absolute failure criterion is shown in Fig. 4.10, and 39.6 ktaps is given as the mode ifi threshold, which is far too early.  Using a 10 equal element model of the beam, the pull-ins and times of failure of the blast loaded beam obtained from the two analyses incorporating different interaction failure criteria are plotted versus the impulses in Figs. 4.11 and 4.12, respectively. As shown in Fig. 4.11, the pull-in versus impulse curves exhibit similarities to the curves of mid-span permanent deformation versus impulse. This is intuitive because the beam will be pulled in more as the beam deforms more severely laterally, provided the beam is inextensible or has limited extensibility. In Fig. 4.12 the first predicted time of failure for both curves occurs around 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 the duration of the pressure pulse).  This also agrees very well with the experimental  observation that mode ifi failures occur early. However, it is clear that the times of failure predicted by the linear interaction criterion are far earlier than the ones predicted by the quadratic interaction failure criterion.  85  0.9 0.8 0.7 0.6  0.4 0.3  -  -  -  -  -  0.2 I  0.1 20  30 _—  I  40 Impulse (ktaps)  tensile strain ratio (10 elem)  -_.._  I  I  50  60  70  shear stress ratio (10 elem)  Fig. 4.10. Comparison of tensile strain ratios and shear stress ratios vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beam using the linear interaction criterion.  86  0.7 0.6 0.5  .  0.4  C  0.3 0.2 0.1 0  0  10  20 ._.  30 40 Impact (ktaps)  QIC (10 elemets)  __  50  60  70  LIC (10 elements)  Fig. 4.11. Comparison of pull—ins vs. impulses obtained from a 0.375 in. by 1 in. by 8 in. beam using the quadratic interaction criterion (QIC) and the linear interaction criterion (UC) with a ten element mesh.  87  350 300  aa) 250 Cl)  2C)  200  4l50  -  100 LIC  50 00  10 _.  20  40 30 Impulse (ktaps)  QIC (10 elements)  ....  50  60  70  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 using the quadratic interaction criterion (QIC) and the linear interaction criterion (LIC) with a ten element mesh.  88 The same failure analysis procedure incorporated with two different interaction failure criteria has also been carried out for a 0.25 in. x 1 in. x 8 in. beam subjected to different impulse loads. The mid-span permanent deformation predicted by the two analyses using a 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 permanent deformation very well. However, the linear interaction failure criterion gives the mode II threshold as early as 17.8 ktaps (from Fig. 4.13) and mode ifi threshold as early as 28.7 ktaps (from Fig. 4.14). According to the experimental data, 28.7 ktaps is reported to be the mode II threshold and 46 ktaps the mode ifi threshold. Comparatively, the quadratic interaction 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 the tensile 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 calculating the response for the case of a 0.375 in. x 1 in. x 8 in. beam subjected to an impulse of 17.8 ktaps, which (from the experiments) is supposed to cause mode I failure.  The time  histories 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 four different meshes, i.e. 8, 10, 16, and 20 equal element grids, are presented. The four time histories are observed to be very similar, and the differences between the curves decreases as the mesh becomes denser. That is, the maximum strain at the support is independent of the mesh and converges as the mesh is refined. Therefore, the tensile influence portion of the failure condition, which is described by tensile strain ratio refining of the mesh.  will converge with  89  2.5  C C 0  2  E I  a)  4-.  C  a)  C Cu  a)  0 C  .O.5  0  0  .10  20  30 40 Impulse (ktaps)  QIC (10 elem) A  Experimental Results  -_-_  x  50  60  70  LIC (10 elem) Shen & Jones’  Fig. 4.13. Midspan permanent deformations vs. impulses obtained from a 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction criterion and the linear interaction criterion with a ten  element mesh.  90  1 0.9 0.8 0.7 0.6  0.4 0.3 0.2  -  -  -  0.1 0 10  20 ——  30  tensile strain ratio (10 elem)  40 Impulse (ktaps) .÷  50  60  70  shear stress ratio (10 elem)  Fig. 4.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 using the linear interaction criterion.  91  1.1 1 0.9 0.8 0  0  0.7 0.6 0.5 0.4 0.3 0.2  20  30  40 Impulse (ktaps)  tensile strain ratio (10 elem)  .  50  60  70  shear stress ratio (10 elem)  Fig. 4.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 using the quadratic interaction criterion.  92  0.09  0.08  0.07 t 0 0.06 Cl)  0 .  0.05  0.04  0.01  0  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  Time (millisec.) —.-—  8 elements  _-_-  10 elements  —i---  16 elements  —--  0.9  1  20 elements  Fig. 4.16. Time histories of total strain at the support for different finite element meshes.  1.1  93  6 5 4 3 C  o  0  o2 cO  z Cu Q)0 h  —2 —3 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  Time (millisec.) ——  8 elements  ——  10 elements  —--  16 elements  ——  20 elements  Fig. 4.17. Time histories of wall reaction for different finite element meshes.  1.1  94  Comparisons of wall reactions calculated from different finite element meshes are presented in Fig. 4.17. The same four meshes as used above are employed. During most of the time history, the four curves are almost the same. At the very early stage of the response (0  -  75 sec) there is a small difference, but as the mesh is refined the results  converge. The results from the 16 element grid are very close to the results calculated from the 20 element grid. When the stable response starts (500 j.i.sec) there is a small shift between the responses calculated from the different meshes. This is probably because the principal frequencies of the different finite element lumped mass models of the beam are slightly different.  Yet, as shown in the figure, the time history of the wall reaction  obtained from 16 element mesh is almost exactly the same as the one from the 20 element mesh. That is, with refining of the mesh, the results converge. From the point of view of the failure criterion, it is the wall reaction force in the early state of response which plays an important role.  As previously discussed, the results of wall reaction force are  independent of the element size and converge when the mesh is refined. Therefore, the shear stress ratio  / tuli  is independent of mesh size, as is the shear portion of the  failure criterion.  Summarizing the above discussion, both tensile and shear parts of the failure criterion are independent of the mesh and converge as the mesh becomes denser, and therefore the prediction 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 force help to explain Fig. 4.9. This figure, during the second failure mode, shows the tensile strain ratio and shear stress ratio as almost the same regardless of the mesh.  This is  because the maximum strain at the support, which is dominant, is independent of element size and the wall reaction force is almost the same for different meshes between 75  -  500  95 sec. When the beam starts to fail in mode Ill, the ratio curves start to be different. This is because although the wall reaction force (the controffing factor in mode ifi failure) is also mesh independent, there are differences in results calculated from different meshes in the 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 ratio converge as the mesh is refined during the mode ifi failure part.  2. Another Two Interaction Failure Criteria As discussed in Chapter 2, from variational principle, the relation between the resultant forces and wall reaction force is obtained. 6W  =  71——  ã?73  +V  16 1 6W 16 —pI’iJ + —pI—-—lCdIW ..  -  .  1  6W  (4.2)  Thus, in a case when the wail reaction force is not available from the equilibrium (such as an 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 an impulse of 17.8 ktaps. In this case, as stress related results are more accurate in middle of the element, the resultant forces at the middle of the boundary element are used to calculate F, the transverse shear resultant force transmitted along the length of the beam, to approximate the wall reaction R.  6x  ax  (4.3)  Because the midpoint of the boundary element is quite close to the boundary, the acceleration and velocity related terms are neglected in the calculation. As the element mesh is denser, the F should be able to approach the wail reaction. Comparisons of F (the transverse shear force) and the wall reaction force are presented in Fig. 4.18 for an 8 element mesh and Fig. 4.19 for a 20 element mesh. As observed from the graphs, for both  96  6 5 4 3 .0 Cl)  w  o2  —2 —3 0  0.1  0.2  0.3  Fw  0.4  0.5  0.6  Time (millisec.) —---  0.7  0.8  0.9  Wall reaction  Fig. 4.18. Comparison of time histories of Fw and wall reaction for an eight element mesh.  1.1  97  5  4-  3-  2 0  0)  D  h  -  —2  I  3  0  0.1  I  0.2  0.3  Fw  0.4  I  0.5  0.6  0.7  0.8  Time (millisec.) __  0.9  1  Wall reaction  Fig. 4.19. Comparison of time histories of Fw and wall reaclion for a  twenty  element mesh.  1.1  98 meshes the Fs at the mid point of the first element are very close to the wall reaction forces. For a denser mesh (e.g. 20 element mesh) the approximation is slightly better. As previously discussed, the wall reaction will converge as the mesh is refined. Hence the resultant 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 along the beam at seven different times. The beam and the loading condition are the same as mentioned above. A 20 element mesh is used to model the beam.  It can be clearly  observed how the shear forces develop during the loading. When the loading stops at 60 p.sec, the shear force at the boundary reaches its maximum.  Using the interaction failure criterion with the shear stress calculated from the finite element analysis, analyzing the 0.3 75 in. x 1 in. x 8 in. beam subjected to the whole range of impulsive loads results in the mid-span permanent deformation versus impulse figure as shown in Fig. 4.21. As before, the linear interaction failure criteria predicts the mode II threshold earlier than the quadratic interaction one, and underestimates the mid-span permanent deformation during mode II failure.  For the quadratic interaction failure  criterion, results from three meshes are compared in the graph. All three curves fit the data well and predict the increasing-decreasing-constant trend of mid-span permanent deformation. The mid-span permanent deformation improves when the mesh is refined to 12 elements from 10 elements. This is because stresses obtained from the program are used in the analysis and a finer mesh is needed to give accurate enough results, as opposed to an analysis based on displacement related quantities. The mode II threshold predicted by this quadratic interaction failure criterion is 33.4 ktaps, compared with 46 ktaps according to the experiments. Mode ifi threshold can be specified from the comparison of the tensile strain ratio and the shear stress ratio, which is shown in Fig. 4.22. As seen from the figure, after a large increase in the shear ratios, they remain quite constant past  99  6  5  4  3 1  C) 0  C  a  I-  0 Cu  1  C,)  0  —1  —2 0  1  10 microsec.  2  3  Distance from the support On.) —  20 microsec.  __  30 microsec.  4  5  40 microsec.  6 5 4 3 .  Q C)  C  o -  CS  0  1  CS Cl)  0 —1 —2 —3 0  1  40 microsec.  3  2  —_  Distance from the support On.) 50 microsec. 60 microsec. —_  4  —a— 70 mlcrosec.  Fig. 4.20. Distributions of shear force along the half beam in the early stage of response.  5  100  3  .  2.5  C 0  E2  0  10  20  30  40  50  60  70  Impulse (ktaps) —.—  _-_  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. beam using the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) and the linear interaction failure criterion with shear stress from the FE analysis (LICF).  101  1.1  0.6  :  -  she 70  Impulse (ktaps) __  _-  tensile ratio (10) shear ratio (10)  __-  —.  tensile ratio (12) shear ratio (12)  ——  —h--  tensile ratio (14) 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 the FE analysis (QICF).  102 61.6 ktaps. Then this is specified as the mode ifi threshold, which is very close to the experimental result of 64 ktaps.  The mesh is then further refined to 14 elements, with the first element in the 10 equal element mesh subdivided into 4 equal length elements and the second element halved. The results in the mid-span permanent deformation are almost the same as from the 12 element mesh. However, a large difference is observed from the comparison of tensile strain ratio and shear stress ratio (Fig. 4.22) because of the boundary layer effect (discussed in Chapter 2). The distribution of the first order shear force at the midpoint of each element along the beam at 60 Lsec, which is the duration of the load, is shown in Fig. 4.23. As is evident from this figure, the shear force at the midpoint of the first element is smaller than that in the middle of the second element. This is because the midpoint of the first element is inside the boundary layer, within which the shear force erroneously decreases from the maximum to zero at the very beginning boundary of the beam. To avoid this, an element greater 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 to various 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.8 and 39.6 ktaps, which compares well with the 28.7 and 46 ktaps reported from experiments.  3. Failure Criterion Using Theory ofPlastic Work Density To assess the accuracy of the failure criterion using sectional plastic work density in modelling the failure of the blast loaded beam, a series of finite element failure analyses incorporating this failure condition have been done for a 0.375 in. x 1 in. x 8 in. beam subjected to different impulsive loading conditions.  The mid-span permanent  103  5  4  3  —.  a) o  0 .  2  C 0  0°’ h.  a)  a,  0  I-i  Cl, 0  —2 0  1  2  3  Distance from the support (in.)  4  Fig. 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.  5  104  3  • 2.5  .-  0 22  0.5  0 30  Impulse (ktaps) .  QICF (10 elem) Shen & Jones’  —..—  —---  LICF (10 elem)  A  Mode I exp. results  QICF (12 elem)  Fig. 4.24. Midspan permanent deformations vs. impulses for a 0.25 in. by 1 in. by 8 in. beam using the quadratic interaction failure criterion with shear stress from the FE analysis (QICF) and the linear interaction failure criterion with shear stress from the FE analysis (LICF).  105  1.1  tensile strain ratio  0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2  -  -  —  -  -  shear stress ratio -  0.1 0 20  30  40  Impulse (ktaps) _.—  —  tensile strain ratio (10 elem) shear stress ratio (10 elem)  —.--  50  60  70  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 FE analysis.  106 deformations 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 elements of equal length, and 12 elements with eight of equal length close to the mid-span and four of half size near the support.  Results are presented in Fig. 4.26, along with  the  experimental data. As seen from the graph, neither curve predicts the trend of mid-span permanent deformation well. Furthermore, there is quite a large difference between the curves obtained from the two different meshes. Analysis based on the 10 element mesh predicts the mode II failure at 25.6 ktaps, while the 12 element mesh analysis predicts 21.8 ktaps. This suggests that analysis using the failure condition of sectional plastic work density is mesh dependent. This is probably because the sectional plastic work density is calculated at each spanwise Gauss point. The coarser the mesh, the fewer the Gauss points, and thus there are fewer sections to be considered in the failure condition, and the more unlikely the beam will fail. To overcome this constraint, the failure criterion is changed into comparing the plastic work done in the first element (which is checked to be always the highest) to the critical value. The results from the analysis using the changed failure condition employing a 10 element mesh are also shown in Fig. 4.26. The results are still unsatisfactory as the trend of the mid-span permanent deformation is not accurately portrayed and the mode II threshold is predicted as early as 25.6 ktaps. This amendment did not improve the prediction probably because the stress state at the boundary of the beam undergoing large plastic deformation is three dimensional, and the one dimensional beam theory cannot accurately describe the beam behaviour near the boundary. The plastic work density calculated near the boundary cannot represent the real value exactly, making it difficult to use the failure condition based on plastic work density theory to predict the rupture of the beam.  Although the failure criterion based on plastic work density does not predict the failure of the beam satisfactorily, a detailed investigation of the distribution of sectional plastic  107  3  .  2.5  C 0  E  2  0)  C a, 1.5 C a,  E a,  0 C a, 0 U,  0.5  0  0  ..  10  20  30  40  50  Impulse (ktaps)  SPWD (10 elem)  SPWD (12 elem)  Shen & Jones’  integrated SPWD  A  60  70  Experimental Results  Fig. 4.26. Midspan permanent deformations vs. impulses for a 0.375 in. by 1 in, by 8 in. beam using failure conditions based on sectional plastic work density (SPWD).  108 work along the beam will help in understanding the distribution of plasticity in the beam when subjected to different impulses.  Shown in Fig. 4.27 are a series of spanwise  distribution of sectional plastic work density for a 0.3 75 in. x 1 in. x 8 in. beam subjected to an impulse of 17.8 ktaps, which caused mode I failure in the experiments. The sectional plastic work density shown is evaluated at the midpoint of each element.  The same  distributions for the same beam subjected to 42.9 ktaps and 64 ktaps (which caused modes II and ifi failure experimentally) are presented in Figs. 4.28 and 4.29. In the figures, 107 and 46 sec are the times of mode II failure predicted by the quadratic interaction failure criterion for the two cases, respectively, while 120 and 66 J.Lsec are the times of mode II failure 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 the highest. The development with time of the plastic work in the half beam for mode I failure can 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 of the (whole) beam. In mode II, rupture of the beam occurs before much plasticity develops in the middle of the beam (shown in Fig. 4.28). In mode ifi the beam fails even earlier and the 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 then ifi, and the plasticity in the beam becomes more localized. This agrees with experimental observations.  Another detailed study of the time histories at different energies gives some idea about the different energetic characters of different failure modes. The time histories of external work, 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 correspond to modes I, II, and ifi failure of the beam, are presented in Figs. 4.30 (a), (b) and (c)  109 400  300 0  -x 0  :200 0 0  a. 0 C 0  a 100 0  0)  0  Distance from the iupport (en.) _.  25 microsec.  ....÷....  50 microsec.  -__  75 miorosec.  100 mioroseo.  800 700 >.  600 0  ¶so  300 C 0  iS e  200  0) 100  Distance from the Lpport (in.) _..  150 miorosec.  ........  175 microsec.  _-_  200 microsec.  1500  to  C  1000  500  0  Distance from the lupport (in.) 225 microsec.  _,_  250 microsec.  ___  275 microseo.  -__  300 mlorosec.  _-__  325 microsec.  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.  110  2500  >%2000  4-  Cl) C  Cu 0  1500 C) Cu Q-iooo  cu C 0  Cu 500  0 0  1  2  3  Distance from the support (in.) —.--  9 microsec.  ——75 microsec.  —--  —--  25 microsec. 100 microsec.  —-_  _-—  50 microsec. 107 microsec.  Fig. 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.  4  111  1500  >‘  0 C  1000  I o 13  500  Cl)  0  Distance from the support (in.) 6 microsec.  _  25 microsec.  __-  46 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.  112 1200  1000  600  C  600  w  400  200  0  a 7  6  5  [ 3  2  0  20  C.  external work 15  kinetic energy 910  5  plastic work 0  0  0.5  Time (millisk.)  1.5  2  Fig. 4.30. Time histories of different energies for half a 0.375 in. by 1 in. by 8 in. beam subjected to impulses of (a) 17.8 ktaps. (b) 42.9 ktaps. (C) 64 ktaps.  113 respectively. As seen from Fig. 4.30 (a), during the mode I failure the majority of the external work is absorbed into the plastic work in the beam; this corresponds to the plasticity 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 work absorbs the majority of the external worlç this corresponds to the plasticity developed in the beam decreasing. Finally, in mode ifi failure (Fig. 4.30 (c)) most of the external work is turned into kinetic energy of the beam. The plastic work in the beam is relatively very small; this corresponds to the beam failing at an early time and plasticity in the beam is localized.  The residual kinetic energy is further plotted vs. impulse for half a broken beam of 0.3 75 in. by 1 in. by 8 in., as shown in Fig. 4.31. The highest velocity of the flying broken beam is 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 be as high as 13975 in./s (355 mIs), which is faster than the sound speed. Research on the blast loaded circular plates [29] also shows the same trend in the residual kinetic energy vs. impulse relation.  4.4.3 Convergence and Energy Conservation  Having 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 time steps, 1, 0.5, and 0.1 sec, have been used in analyzing the beam of 0.375 in. x 1 in. x 8  114  20  15  -  C .  LIC 0)  5 1O Co W E 1 0 C 5-  QIC I  0 0  10  20  30  I  I  40  50  Impulse (ktaps)  QIC (10 elements)  .  60  70  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 different impulses.  115  30  25  -  -  C  20  LIC  —  2 •0  > J15 C  w 0 C  10  -  5-  QIC 0— 0  10  20  30  40  Impulse (ktaps)  QIC (10 elements)  .  50  60  70  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 different impulses.  116 in. subjected to an impulse of 17.8 ktaps. A 10 equal element mesh is used. The mid-span permanent deformations obtained are listed as follows for comparison:  Table 4.1 Convergence of mid-span permanent deformations with time step. Time step (.isec)  max Aw  mm Aw  Aw  m  m  m  1 0.5 0.1  0.64495 0.64244 0.64034  0.566235 0.55977 0.55758  0.60235 0.60108 0.59896  The difference between results from calculations using 1 and 0.1 sec as time steps is only 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 half beam, 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 of mid-span permanent deformations with mesh. No. of elements 8 10 16  max Aw in.  mm Aw  m  m  0.63716 0.64495 0.64772  0.55488 0.56235 0.56528  0.59602 0.60365 0.6065  From the table it can be observed that the results are converging well.  117 For the energy conservation check, two cases of a 0.375 in. x 1 in. x 8 in. beam subjected to 17.8 ktaps and 52.9 ktaps are analyzed using different time steps. The results are presented in Figs. 4.33 and 4.34 respectively. In these figures the time histories of the energy difference ratio is shown. Energy difference ratio rE is defined as  rE  =  Et -Ek pl Eext  (4.4)  where Ee,, Ek, and E 1 correspond to external energy, kinetic energy, and plastic work done in the half beam. Since the Eext should be greater than the summation of Ek and with 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. As observed from the figure, when rE is less than 0, it is always greater than -1.5%, which is still acceptable.  With the decrease in the time step, the rE is positive, so the results  converge to the correct value. Therefore, the energy is conserved when a small enough time step is used.  4.5 Summary  In this Chapter, after the series of experiments conducted by Menkes and Opat are  introduced, the numerical simulation of the experiments by using finite element analysis incorporated with the failure conditions are presented. The results obtained by employing different failure criteria are compared to the data from experiments in the discussion. The convergence is also presented and, as shown, the results are converging as the time step and mesh are refined.  118  0.1 0.09  -  0.08  o  -  0.07  -  4-  0.06C  0.05  -  a)  0.04> 0.03C  0.02  -  0.01  -  0 I  I  0.4  0.5  —0.01 0  0.1  0.2  0.3  0.6  Time (millisec.) —.--  ...  ts=1 microsec. (10 elem)  __  0.7  0.8  0.9  1  1.1  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 1 in. by 8 in. beam subjected to an impulse of 17.8 ktaps.  119  0.02  0.015  .2  0.01  a) 0.005  a) V  a)  C W —0.005  —0.01  —0.015 0  —i—  0.5  1  Time (millisec.)  1.5  2  ts=1 microsec. (10 elem)  ts=0.5 microsec. (10 elem)  ts=O.1 microsec. (10 elem)  ts=0.1 microsec. (12 elem)  2.5  Fig. 4.34. Time histories of energy difference ratio using different time steps for a 0.375 in. by 1 in. by 8 in. beam subjected to an impulse of 52.9 ktaps. The rupture of the beam is predicted by the quadratic interaction criterion with shear stress from the FE analysis.  120 According to the comparison, the quadratic interaction failure criteria using shear stress from the finite element analysis provided the closest answer, but at the cost of using a denser mesh.  However an extremely dense mesh should be avoided because of the  boundary layer effect of the higher order beam theory solution. The quadratic interaction failure condition using wall reaction force provides a good prediction as well while allowing a comparatively coarse mesh to be used to obtain accurate results; thus this technique is more practical for engineering applications.  From the discussion, the effect of interaction cannot be neglected, particularly in the mode ifi failure. The membrane effect still plays an important role in the shear dominant failure of the beam.  Additionally, the ruptures of the beam predicted from the analysis are earlier than the experimental results, which is desirable for engineering practice.  This conservatism is  probably because the beam ends in the experiments may not be fi.illy restrained according to the configuration of the experiment and the length of plastic hinge cdi might be smaller than in reality.  121  CHAPTER 5 CONCLUSION  5.1 Summary  A failure analysis procedure to predict the rupture of the impulsively loaded ductile clamped beams which exhibit geometric and material non-linearities has been presented in this 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 displacement  relation by using the principle of virtual work and Total Lagrangian Approach, the governing variational equation of motion is obtained. It is yet to be noted that the beam kinematics is restricted to the range of small strains and moderately large displacements. After the finite element spanwise discretization is applied, the variational governing equation of motion reduces to a system of equations of motion which are differential equations with respect to time. The constitutive model which the system of equations of motion are based on considers elastic-plastic isotropic strain hardening material behaviour by von Mises yield criterion and associate flow rule and strain rate sensitivity by Cowper Symonds relationship. Then the differential equations are integrated over the time domain by the central difference method. Lumped mass matrix is used to improve the accuracy and efficiency. A conservative critical time step is determined for stability. The governing  122 differential equations of motion now further reduce to a system of algebraic equations which can be solved easily.  Based on finite element method calculations interaction failure criteria and the failure condition using plastic work density theory are proposed.  The premise behind the interaction failure criteria is to consider the contribution of both tensile tearing and transverse shearing effects to the failure of the beam. Thus it can explain the overlap of mode II and ifi failures between the thresholds.  The beam is  modelled as an elastic-plastic beam with a plastic hinge at the supports. The plastic hinge assumption is employed in calculating the tensile strain. The shear stress can be obtained either from wall reaction forces (from equilibrium) or first order shear force (from the finite element analysis). Tensile strain ratio and shear stress ratio are used to calibrate the influence of the tearing and shearing effect, respectively. Two interaction functions have been investigated, namely quadratic interaction and linear interaction.  The concept on which the failure criterion using the plastic work density is based is that a crack initiates at a point where the density of plastic work reaches a critical value. It is assumed that the beam does not fail until the whole section fails since the rupture always occurs in a section. Therefore the failure conditions can be put that the beam will not fail until the sectional plastic work density reaches a critical value.  Because when the rupture occurs there still remains quite a lot of residual energy, the post failure analysis on a blown off free beam is presented. In this way, the mid-span permanent deformation from the analysis is comparable to the experimental results.  123 To assess the accuracy of the failure conditions, analysis using different filure  conditions have been carried out to simulate the experiments conducted by Menices and Opat. Results obtained from these analyses have been compared to the experimental data. For the quadratic interaction filure conditions, encouraging results have been obtained. It is to be noted that Menkes and Opat admitted in their paper that the specification of mode II and ifi thresholds are highly subjective. Therefore, the conservative results obtained from the analysis are desirable in the engineering sense.  As stated in the introduction, the objective of the present study is to determine a simple failure model to be incorporated into the finite element analysis so that the onset of the mode II and ifi failures of impulsively loaded beam can be predicted. satisfactorily accomplished.  This has been  Furthermore, some insight into the characteristics of the  different failure modes is obtained from the discussion about the numerical results.  5.2 Areas Requiring Further Investigation  In 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 for further investigation.  During the formulation of the finite element analysis of the problem, it is assumed that  the strain is everywhere small. However, when the rupture of the beam occurs due to excessive strain, the strain becomes significant. validity of simplification in the formulation.  This might have some effect on the  124 In the present work, the rupture strain of the beam is chosen to be equal to the static uniaxial rupture strain. The ultimate shear strength is determined based on the von Mises yield criterion and the ultimate tensile strength of the material. That the dynamic shear strength 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 is based on several assumptions which are not confirmed. An investigation is therefore needed to check the sensitivity of the results to these two parameters.  During the modelling of the impulsive loading condition, a rectangular pressure pulse is  used and an assumption was made as to the duration of the load. Since the assumptions on which the load modelling is based are not well confirmed, this introduces the possibility that the modelling of the load shape and duration might have some influence on the failure of the beam. A parameter study is needed to confirm or improve the modelling of the blast loading condition in the future.  While calculating the bending strain at the support, the length of the plastic hinge is defined as l=ah, where cx, the plastic hinge parameter, is chosen to be two. However, as the beam goes into the membrane stage, the length of the plastic hinge will be greater, thus decreasing 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, the length of the plastic hinge will decrease, as indicated in the figure on distribution of plastic work along the beam. This will have an effect on the time of rupture of the beam and thus the mid-span permanent deformation. As discussed above, the specification of the length of the plastic hinge, hence the plastic hinge parameter a, will have an effect on the accuracy of the interaction failure model; possibly refinements are necessary.  125  References  1.  Jones, N. “A Literature Review of the Dynamic Plastic Response of Structures,” Shock and Vibration Digest, Vol. 7, August 1975, . 105 89 pp.  2.  Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 1 and 2, Shock and Vibration Digest, Vol. 10, September 1978, pp.21-33 and October 1978, . 19 13 pp.  3.  Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 3, Shock and Vibration Digest, Vol. 13, October 1981, . 16 3 pp.  4.  Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 4, Shock and Vibration Digest, Vol. 17, February 1985, . 47 35 pp.  5.  Jones, N. “Recent Progress in the Dynamic Plastic Behaviour of Structures”, Part 5, Shockand Vibration Digest, Vol. 21, 1989, . 13 3 pp.  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 of Engineeringfor 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, . 159 133 pp.  9.  FoIz, B. R., “Numerical Simulation of the Non-Linear Transient Response of Slender Beams”, MA.Sc. thesis, University of British Columbia, Vancouver, Canada, April 1986.  126 10. Fagnan, J. K., “Failure Prediction of Blast Loaded Ductile Beams Using Approximate Deflection Profiles”, Ph.D. Thesis Proposal, University of British Columbia, Vancouver, Canada, June 1991. 11.  Shen, W. Q. and Jones, N., “A Failure Criterion For Beams Under Impulsive Loading”, 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 Metal Beams Subjected to Large Dynamic Loads”, Structural Crashworthiness and Failure, N. Jones and T. Wierzbicki, eds., Elsevier Applied Science, 1993, pp.95130. 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, Pergamon Press, 1982. 15. Levinson, M., “A New Rectangular Beam Theory”, Journal of Sound and Vibration, Vol. 74, 1981, pp.81-87. 16. Bickford, W. B., “A consistent Higher Order Beam Theory”, Developments in Theoretical andAppliedMechanics, Vol. 11, 1982, . 150 137 pp. 17. Heyliger, P. R. and Reddy, 3. N., “A Higher Order Beam Finite Element for Bending and 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, . 752 745 pp. 19. Reddy, J. N., Energy and Variational Methods in Applied Mechanics. New York: John Wiley. 1984. 20. Palmer, A. C., Structural Mechanics, Clarendon Press, 1976.  127 21. Cook, R. D., Malkus, D. S. and Plesha, M. E., Concepts and Applications ofFinite Element 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 Element Method, University ofBritish Columbia, Vancouver, BC, January, 1993.  24. Olson, M. D., CIVL 537 Course Notes: The Finite Element Method, University of British Columbia, Vancouver, BC, September, 1992. 25. Jiang, J., and Olson, M. D., Documentation for FENTAB Department of Civil Engineering, University ofBritish Columbia, Vancouver, 1990. 26. Nonaka, T., “Some Interaction Effects in a Problem of Plastic Beam Dynamics, Parts 1-3”, Journal ofAppliedMechanics, Vol. 34, 1967, . 643 623 pp. 27. Fagnan, J. R., FENTABF, Department of Civil Engineering, University of British Columbia, 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 and Tearing ofBlast Loaded Circular Plates”, to be published.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050411/manifest

Comment

Related Items