A Review of Soft Body Impact on Composite Structure

Share

Transcript

1 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference18th AIAA 2010-2865 12 - 15 April 2010, Orlando, Florida A REVIEW OF SOFT BODY IMPACT ON COMPOSITE STRUCTURE Arafat I. Khan1 and Rakesh K. Kapania2 Virginia Polytechnic Institute and State University, Blacksburg, VA, 24061, USA Eric R. Johnson3 Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA Recent progress in material modeling and numerical analysis related to soft body impact on fiber reinforced composite structure is presented. Material modeling of the composite includes continuum damage mechanics for intra-ply progressive failure and cohesive zone models for interface fracture. Soft body projectiles exhibit significant deformation on impact and flow over the structure spreading the load. A numerical model of the soft body that has been correlated to the test data is the smooth particle hydrodynamics method in conjunction with an equation of state. The composite material model and the soft body projectile model are implemented in explicit finite element software packages, which also contain contact algorithms and non-linear transient analysis capabilities. Nomenclature = strain vectors = stress vectors S = elastic compliance matrix E1, E1, E12 = Youngs moduli for composite ply G12 = shear modulus for composite ply d1, d2, d12 = damage parameters for the composite degradation model Y1, Y2, Y12 = damage functions for composite degradation f1, f2, f12 = damage evolution functions e = elastic strain p = plastic strain = effective shear stress R = isotropic hardening function p = accumulated effective plastic strain = yield stress 33 = tensile stress u3 = displacement jump across the interface in cohesive zone model k3 = the tensile stiffness u30, u3m = displacements at the peak stress 33m = maximum peak stress for delamination model GIC = critical fracture energy C0, C1 = material constants C2,C2 = material constants = ratio of current density and initial density = current density 0 = initial density vs = shock velocity 1 Graduate Research Assistant, Aerospace and Ocean Engineering, [email protected] 2 Mitchell Professor, Aerospace and Ocean Engineering, [email protected] , Associate Fellow, AIAA 3 Professor Emeritus, Aerospace and Ocean Engineering, [email protected], Senior Member, AIAA 1 American Institute of Aeronautics and Astronautics Copyright 2010 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.

2 vp = particle velocity c0 = speed of sound for water i = center node for SPH method j = neighboring nodes in SPH method rij = distance between nodes di, dj = diameter of the center node and neighboring node Mi , Mj = mass of the center node and neighboring node = smoothing distance for SPH method = smoothing function = volumetric strain = normal strain rates = shear strain rate = volumetric strain rate = normalization parameters for smoothing function = stretching radial velocity = nodal force in the x direction = force due to the in plane stress = force due to the hoop stress = force in z direction 0 = the maximum allowable overlap distance of an SPH node and a slave node I. Introduction C OMPOSITE materials are increasingly used in many aircraft structures due to their high specific strength, light weight, resistance to fatigue/corrosion and flexibility in design. Research based on a soft body impact on a composite structure by a bird, hailstones and tire rubber are of great interest in the aerospace industry. A soft body impact on an aircraft structure such as the fuselage or the engine can lead to serious damages to the aircraft. Bird strike contributes significant economical damage in aerospace industry. According to Canadian Transport[1], The problem goes back to the early stages of aviation when the Wright brothers recorded the first bird strikes on September 7th 1905 It is estimated that a bird strike event occurs once in every 2000 flights. According to Georgiadis et.al. (2008), in civil aviation more than 50 planes and 223 lives are reported to have been lost due to bird strikes since 1912 and the cost of damage due to bird strike in worldwide aviation industry is about 1 billion USD per year. The most recent event of birds striking aircraft engines occurred on January 15 th, 2009, when a US airways aircraft collided with a flock of birds. Both engines lost power as a result of the impact, and the aircraft was forced to land into the Hudson River in New York City. An important monograph by Abrate (1998) presents research on the impact of composite structures. Topics in this monograph include how damage develops during impact, the effect of impact-induced damage on the mechanical behavior of structures, and methods of damage prediction and detection. However, the specific topic of soft body impact on composites is not discussed. Models of the soft body projectile are, of course, much different than models of a hard body projectile. Soft bodies have significant deformation on impact and flow over the structure spreading the impact load. A bird is usually represented using a substitute material, e.g. gelatin. Generally, soft body impact testing is conducted by shooting a gelatin projectile from a gas gun at several different velocities at a target composite structure, Georgiadis et al. (2008) and Roberts et al.(2005). The targets in these tests were composite plates, half rings, and a movable trailing edge. In the following sections of this review we present progressive failure models of the composite material, the smooth particle hydrodynamics model of the soft body projectile, experiments, and commercial finite element software used in modeling soft body impact on composite structures. II. COMPOSITE MATERIAL DEGRADATION MODEL The continuum damage mechanics (CDM) failure models as studied by Talreja (1994) are described in this section. Several composite degradation model was put together by Soden et al.(2004) in the world wide failure exercise. In the work of Puck and Schurmann (2003) depicts failure criteria for composite structure mostly based on the experimental approach. The Puck and Schurmann failure criteria was derived from the quadratic failure criteria proposed by Hashin(1980). The amount and the type of failure mechanisms activated depend on some factors like, mass and velocity of the soft body projectile (impact energy level), geometry of the projectile, geometry of the 2 American Institute of Aeronautics and Astronautics

3 structure, type of fiber and/or matrix used for manufacturing of the composite structure and stacking sequence of the plies. Different kinds of failure modes have been used to model impact on a composite structure, as found in the work of Johnson and Holzapfel (2003). A. Elastic Ply Damage Model The fabric reinforced composite ply is modeled as a homogeneous orthotropic elastic, or elastic- plastic, material whose properties are degraded on loading by microcracking prior to ultimate failure. A continuum damage model formulation employs ply degradation parameters as internal state variables which are governed by damage evolution equations. Constitutive laws for orthotropic elastic materials with internal damage parameters take the general form: (1) where and are the vectors of stress and strains, respectively, and S is the elastic compliance matrix. The elastic compliance matrix is defined in Eq. (2) (2) where d is the damage parameter. The in-plane stress and strain components are, (3) (4) The ply model introduces three scalar damage parameters d1, d2, and d12, where 0

4 (8) where it is assumed that R(0) = 0 and that is the initial threshold value for inelastic stress behavior. The condition F < 0 corresponds to a state of stress inside the elastic domain where the material may be damaged due to elasticity. So it can be deduced that the plastic strain p, (9) showing that in the mathematical model p is the accumulated eective plastic strain over the complete loading cycle. The model is completed by specifying the hardening function R(p). This is determined from cyclic loading tests in which both the elastic and irreversible plastic strains are measured. The general form of the hardening function becomes: (10) The shear plasticity model depends on the parameters , the power index m and the yield stress . C. Delamination Model Delamination failures occur in composite structures under impact loads due to local contact forces in critical regions of load introduction and at free edges. Weak resistance to delamination is a result of the low, resin dominated, through-thickness shear and tensile properties found in laminated structures. Cohesive zone models are typically used to represent delamination initiation, growth and arrest. A cohesive zone model of the interface initially has zero thickness, and across the interface there is continuity of interfacial tractions but jumps in displacements. The equations of the model are given here for the case of mode I tensile failure at an interface. Let 33 be the tensile stress applied at the interface, u3 the displacement jump across the interface and k3 the tensile stiffness. The cohesive zone stressdisplacement model is assumed as (11) For (12) in which the tensile damage parameter is denoted by d3, and c3 = u3m/ (u3m - u30). It can be verified that with this particular choice of damage function d3, the stressdisplacement function has the triangular form shown in Figure 1. Displacements u30 and u3m correspond to the displacements at the peak stress 33m and at ultimate failure respectively. The damage evolution constants are defined in terms of 33m and GIC, the critical fracture energy under mode I interface fracture, by u30 = 33m /k3 and u3m = 2GIC / 33m. Figure 1: Idealized mode I interface stress-displacement function From these expressions, it can be shown that the area under the curve in Fig.1 above is equal to the fracture energy GIC. This interface model therefore represents an initially elastic interface, which is progressively degraded after reaching a maximum tensile failure stress 33m so that the mode I fracture energy is fully absorbed at separation. The process is irreversible for 0 < d3

5 GI = (13) The strain energy release rate is monitored and, if this is found to exceed the critical fracture energy value then the crack is advanced. For mode II interface shear fracture a similar cohesive zone model to mode I is assumed, but with equivalent set of damage constants, and critical fracture energy GII. In general, there will be some form of mixed mode delamination fracture involving both shear and tensile stresses. This is incorporated in the model by assuming a mixed mode fracture condition, which for mode I/mode II coupling could be represented by an interface failure envelope such as Crisfield et. al (1997): (14) where GI and GII are the monitored interface strain energy release rates in modes I and II respectively, GIC and GIIC are the corresponding critical fracture energies, and the constant n is chosen to t the mixed mode fracture test data. Typically n is found to between 1 and 2. Failure at the interface is imposed by degrading stresses when, 1 and the corresponding shear relation. When, 1 there is delamination and the interface separates. D. Fiber Metal Laminates McCartly, et al. (2004) model the leading edge of a wing containing metal plies as well as fiber/matirx plies with respect to bird impact. In order to perform the analysis, a failure model for the metal layer of the composite is required as referred by McCarthy et al. (2004). The material law used for the metal layers is an isotropic elastic plastic thin shell model, with a plasticity algorithm that includes transverse shear effects. In addition, an isotropic damage law has been added to the basic elasticplastic formulation to allow cut-off of element strength at high strains. The elasticplastic behavior was input as a power law: (15) where a is the yield stress, p is the effective plastic strain, and b and n defines the shape of the curve. An isotropic damage law is incorporated in the analysis which states that the total element stresses at each integration point are acted upon in the following fashion: (16) where is the damaged full stress tensor, d is the isotropic scalar damage function, is the plastic strain and is the full stress tensor as calculated from the undamaged elasticplastic material law. III. Modeling the Soft Body Projectile A. Physical Modeling of Projectiles In order to perform analysis of soft body impact on composite structure, the physical modeling of a soft body weighing 4 Lbs (1.82 Kg) is required. The physical characteristic of the soft body results in different impact characteristic. At low velocity, the characteristic of the soft body like a bird affects the impact characteristics. The bird tissue plays an important role at low velocity impact. For high velocity impact cases, these characteristics of the soft body can be neglected. The geometric configuration of the soft body has a significant effect on the impact with a composite structure. Three different geometric configurations for modeling a soft body is discussed by Meguid et al (2007). These geometric configurations are the straight-ended cylinder, hemispherical-ended cylinder, and the ellipsoid. The three different geometric configurations for modeling the soft body projectile are shown in Figure 2. 5 American Institute of Aeronautics and Astronautics

6 Figure 2: (a) straight-ended cylinder, (b) hemispherical-ended cylinder, and (c) ellipsoid. The geometric configuration of the impact projectile has a strong influence on the contact area developed between the projectile and target. The parameters L and D shown in the Figure 2 determine different aspect ratio for the soft body. The geometric configuration and the aspect ratio affect the pressure distribution and impulse response of the impact. B. Smooth Particle Hydrodynamics (SPH) Model In the following paragraphs, we summarize the smooth particle hydrodynamics (SPH) method to model the soft body projectile during impact. More details of the SPH method is provided in the Appendix. The impact between any soft body and structure is considered as a fluids- structure interaction problem. Soft bodies are mostly modeled using gelatin as a substitute bird to analyze soft body impact with a composite structure. Soft body impact limits the local impact damage as the soft body spreads over a significant surface area. Mathematical representation of the material behavior and the structural response for the soft body is a hurdle to modeling a soft body. In the work of Johnson and Holzapfel (2003) it is found that Smooth Particle Hydrodynamics method was used to model the soft body projectile along the equation of state (EOS). Smooth Particle Hydrodynamics is a grid-less computational method whose foundations are in interpolation theory. The material is represented as a set of discrete particles, or interpolation points, which are topologically independent from each other. The two main features of the model are the kernel approximation and the particle approximation. The continuum mechanics conservation equations of mass, momentum and energy are transformed from partial deferential equations to integral equations through the use of a smoothing kernel function which provides the Kernel estimate of the field variables at a position in space. These smoothing kernels must vanish, or become very small outside a radius which is proportional to the smoothing length h, and in the limit that h 0 the kernels approach the delta function. Thus the kernel function defines a range of influence of a point in the continuum. Integration by parts of the integral equations allowed the unknown spatial derivatives of the field variables to be replaced by known derivatives of the kernel functions. The particle approximation is now applied to replace the continuous domain of influence at a point by a set of discrete neighboring particles. It then follows that the transformed integral field equations are replaced by summations over discrete particles within the region of influence of a point. Each particle has an associated mass, velocity and stress-state which will evolve according to the discretized field equations. The development of the discretized field equations is given in Monaghan (1992) along with discussions of suitable kernel functions. The equations of SPH method have been incorporated as the SPH solve option in PAM-SHOCK software package. The SPH solver was incorporated into the PAM-SHOCK in early 1990s as reported by Groenenboom (1997). The SPH solve was incorporated in PAM-SHOCK in order to predict the pressure peak pressures behind a 1-D shock front. C. Equation of State (EOS) in the SPH Model In order to use the SPH model for soft body impact on a composite structure, a constitutive law for the soft body is required with suitable number of materials parameters. The material model used is referred to as an elastic plastic hydrodynamic solid, which was originally developed for ballistic impact in metals, and describes an isotropic elasticplastic material at low pressures, with an EOS describing the hydrodynamic pressurevolume behavior at high pressures. To model soft body impact, the elasticplastic contribution to the material of the soft body behavior may be neglected so that the model reduces to the EOS for the pressure p, which is assumed in this case to have the polynomial form: (17) In these latter equations C0, C1, C2 and C3 are materials constants and is a dimensionless parameter defined in terms of the ratio of current density to initial density 0. This polynomial form is an established approximation of 6 American Institute of Aeronautics and Astronautics

7 the observed EOS for many materials, which is found in the work of Hermainer and Thomas (1995). The polynomial form reduces to a dilatational elastic material law with bulk modulus C1, when C0= C2 = C3= 0. The material constants represent the dynamic behavior of the soft body at impact pressures, and they are difficult to measure directly. Hence, the material constants have to be determined indirectly. Material parameter determination for the polynomial EOS is addressed in the work of Johnson and Holzapfel (2003). The approach used is to calibrate the material parameters by correlating an impact simulation with impact pressure pulse measured in testing. The work of Wilbeck (1977) discusses impact pressures of several soft body soft body projectiles, including rubber, gelatin and chickens, and found that pressure pulses have a certain characteristic form. This form consists of a high peak pressure caused by shock wave propagation in the soft body projectile, followed by a lower fairly constant pressure due to steady flow of the projectile over the target. Furthermore, it was found that the EOS for water could be used as a basis for predicting peak pressures behind a 1-D shock front which were similar to those measured in the test program. For a material such as water which exhibits the linear Hugoniot relation between shock velocity vs and particle velocity vp : (18) where k is a material constant and c0 the speed of sound in the material. Wilbeck shows that the pressuredensity relation across a shock has the general form (19) A series expansion of the latter equation followed by a comparison to the polynomial EOS, identifies the material constants in the polynomial form as (20) (21) (22) The material constant C0 is related to the initial equilibrium pressure in the soft body projectile, and is generally assumed to be negligible. Material constant C1 is bulk modulus of the projectile material. Wilbeck showed that using the values for the Hugoniot constant k and speed of sound c0 for water, that the predicted pressure peak across the shock front for gelatin projectiles fitted test data fairly well. However, for real chickens the pressure pulse had lower peak pressures and was wider. Since the aim of the simulation tools is to simulate real bird impacts on structures, it was recommended to use materials constants in the EOS which represent a mixture of water with about 10 % air. The air content has the effect of reducing density, and lowering the bulk modulus and speed of sound. Thus it is possible to determine the EOS constants using a rule of mixtures model, as proposed by Wilbeck (1977) and Iannucci (1998). This was the procedure used in order to calibrate the SPH model by finding suitable values of p0, c0 and k to represent the EOS of a water/air mixture which fits available impact test data. IV. Experimental Approaches Different experimental approaches to study soft body impact on composite structures are found in the works of several authors. A standard experimental approach includes using a soft body projectile weighing 1.82 kg (4 lbs) as the soft body, Roberts et al. (2005), Georgiadis et al. (2008) and Johnson and Hozapfel (2003). The soft body in the experiments is gelatin wrapped in nylon bag, or gel-pack. According to Georgiadis et al. (2008) three different speed ranges for the gel-pack are used to target the stationary structure. For the lower speed ranges the impact does not result in any significant damage of the composite structure and the gelatin-nylon bag disintegrates. For the medium speed range, the gel-pack produces significant damage in the structure. The high speed test resulted in the soft boy puncturing through the composite structure. The composite structure experiences significant amount of damage due to the soft body projectile impact. Despite the complexity of the test (primary and secondary impacts), the simulation provides excellent agreement to the test. The maximum speed of the gel-pack is used as 225 m/s, and at this impact speed both gel-pack and the composite structure experience large, visible distortions. 7 American Institute of Aeronautics and Astronautics

8 V. A Design to Increase Impact Resistance The development and validation for the design of an innovative composite leading edge structure is described by Kermandis et al. (2005), and is reviewed in this section. The design concept is based on the absorption of the major portion of the bird kinetic energy by the composite skins, in order to protect the ribs and the inner leading edge structure from damage, and thus preserving the tail plane functionality for safe landing. The leading edge skin is fabricated from specially designed composite panels, called tensor skin panels. Tensor skin panels are comprised of folded layers, which unfold under the impact load and increase the energy absorption capability of the structure. It is assumed that the impact load on the composite structure due to a bird should not exceed more than 245 kips. The first approach for solving the soft body impact problem was to study the geometric configuration and the physical parameters involved in the tensor skin design. The tensor skin is a layer of plies within the plies of the skin of a composite structure. A composite skin containing tensor plies has three layers of plies. These plies are called carrying plies, tensor plies and cover plies. These three kinds of plies including the tensor plies are shown in Fig. 3. Figure 3: Stack of plies containing the tensor plies From Kermandis et al. (2005), (with permission) The carrying plies transfer the normal operational aerodynamic load to the ribs and they provide the usual structural skin stiffness. The second group of plies, the tensor plies, contains folded loops between the ribs, which are to unfold when a relatively high lateral load is applied to the skin, e.g. in the case of bird strike. Under normal flight conditions the tensor plies are without any function. In case of a relatively high lateral load all plies except the tensor plies fail. Under the induced penetration load, the tensor plies start unfolding, acting as plastic hinges, until they are loaded in tension like a membrane. Fibers of the composite fabric material unfold from their loop and in a sense catch the projectile. When the tensor plies are completely unfolded, the vertical load component is transferred to the ribs as compressive load. The third group of plies has the purpose to cover the tensor skin. The skin to rib attachment is realized with fasteners conforming to common aviation types of attachment. It is apparent that the functionality of the tensor skin concept as energy absorbing element is based on the capability of the tensor skin material to unfold and therefore high strain-to-failure capability is required. Ordinary composite material systems are too brittle to fulfill this high strain-to-failure requirement. Therefore special material systems had to be considered. A special polyethylene fiber called Dyneema was chosen for the tensor skin concept. Dyneema is a High Performance Polyethylene fiber (HPPE), with high specific strength and modulus as discussed in Ubles et al.(2003). Dyneema fiber fabrics are impregnated by epoxy resin and the resulting composite system attains strains to failure of approximately 3-5% in tension-compression, and of approximately 30% in shear. These failure strain magnitudes make Dyneema suitable for use as the tensor skin material. During production, polyethylene with an ultra high molecular weight (UHMW-PE) is used as the starting material. In the production process of Dyneema, the molecules are disentangled to create filaments; the fiber is drawn and a high level of macromolecular orientation is acquired. This results in a fiber with high specific strength and modulus. Dyneema fiber fabrics are suitable for impregnation by epoxy resin and for laminating within a conventional composite laminate. The limited adhesion of the polymer ethylene fibers to the epoxy contribute to the successful unfolding mechanism of the tensor skin concept. The design process for the tensor skin initially includes determining the energy absorbed by the skin due to the bird impact. In order to apply the tensor skin concept various tensor loop configurations were considered. The tensor load level depends on the tensor angle and therefore on unfolding depth and rib pitch. The total bird kinetic energy has to be absorbed within the leading edge volume; there are limitations to the unfolding depth, which is a function of the tensor skin elongation and the rib pitch. The maximum unfolding depth must fit within the available leading edge space determined by the predetermined leading edge contour. The kinetic energy is determined using, where m is the mass of the soft body and v is the impact velocity. 8 American Institute of Aeronautics and Astronautics

9 Different tensor strips with different loop configurations were considered. A specimen of the tensor strip is shown in Fig4. Figure 4: A small scale tensor strip from Kermandis et al. (2005), (with permission) The simple tensor loop specimens are compared regarding the failure load of the tensor layer, the failure unfolding depth and the specific strength of the ply. Geometric parameters of a tensor loop are show in figure 5. Figure 5: Tensor loop geometrical parameters from Kermandis et al. (2005), (with permission) The larger unfolding depth indicates a more effective energy absorbing concept. The failure load level gives information about the energy that has been absorbed by the tensor layer during a soft body impact. The various kinds of tensor loop configurations considered are listed below in Fig 6. Figure 6: Various tensor loop configurations from Kermandis et al. (2005), (with permission) Small scale tests were performed on the different configurations shown in Fig. 6 to finalize the proper configuration for the soft body impact resistance. According to Kermandis et al. (2005) when a soft body such as a bird strikes the composite tensor skin, the cover ply will fail but the tensor skin will not fail as shown in Fig. 7. Figure 7: Functionality of a Tensor Loop during Impact from Kermandis et al. (2005), (with permission) After the impact the tensor skin will be unfolded and will perform the functionality of a skin to the composite structure. 9 American Institute of Aeronautics and Astronautics

10 VI. Commercial Finite Element Software In order to simulate the soft body impact on a composite structure explicit finite element software packages such as PAM/CRASH/SHOCK [16], LS-Dyna, ABAQUS/Explicit,and MSC.Dytran are appropriate to consider. A review of finite element software for composite impact analysis is presented by Nguyen et al. (2005). The soft body projectile in PAM/CRASH is modeled by Smooth Particle Hydrodynamics method (SPH) along with the appropriate equation of state, Groenenboom (1997). The explicit finite element method is efficient for short duration impact events. Implicit finite element methods are best suited for dynamic analysis where linear or moderate nonlinear effects are to be investigated. The advantage of the explicit FE method is that it easily incorporates geometric and material nonlinearities since the nodal force vector is evaluated at the previous time step to estimate the response at the current time step. Very small time steps are required to achieve an accurate numerical solution. Fewer computational resources are needed, but the very small time step means that only short duration events can be analyzed in an acceptable length of time. Most of the papers reviewed used the software PAM/CRASH and PAM/SHOCK for the numerical simulation of soft body impact on composite structure. PAM/CRASH is impact simulation software developed by the ESI group of France [19]. PAM/SHOCK has been used by many authors to model the projectile for a soft body impact. PAM/SHOCK is a nonlinear Finite Element software package which is used to simulate shock transmission and subsequent transmission responses within materials. PAM/SHOCK addresses laws of shock waves propagation within various environments, it can implement Smooth Particle Hydrodynamics (SPH) and PAM-SHOCK offers a wide range of materials to calibrate behavior under various impact conditions such as fusion, composites deformation etc . VII. Concluding Remarks This review presents important features required for an accurate numerical simulation of a soft body projectile impacting a composite material structural component. Laminated composite material degradation needs to include continuum damage models for intra-ply failure and cohesive zone models for interfacial fracture initiation and growth. Soft bodies have significant deformation on impact and flow over the structure spreading the impact load. Hence, the soft body projectile is modeled using the smooth particle hydrodynamics method in conjunction with an equation of state. The projectile and composite material models are implemented in explicit finite element software packages such as PAM/CRASH, LS-Dyna, MSC Dytran, and ABAQUS/Explicit. The commercial software package also needs to contain contact algorithms and the capability to analyze nonlinear transient response. Tests to validate the numerical simulation of soft body impact on a composite are provided in the papers by Georgiadis, et al. (2008) and Johnson, et al. (2003). In addition, we cite the paper by Kermandis et al. (2005), which describes the design of an innovative composite leading edge structure to absorb the impact energy by unfolding composite plies to catch the projectile. Appendix Smooth Particle Hydrodynamics SPH (Smooth Particle Hydrodynamics) techniques provide the capability to perform high distortion impact computations in a Lagrangian framework as mentioned in Johnson et al. (1996) and Monaghan (1992). The basic SPH technique was first introduced by Lucy (1977) and Gingold and Monaghan (1977),and two comprehensive reviews are presented by Benz (1989) and Monaghan (1992).The effect of strength was added by Libersky and Petschek (1990), axisymmetric algorithms were developed by Johnson et al.(1993) and Petschek and Libersky (1993) and a Normalized Smoothing Function algorithm was developed by Johnson and Beissel (1996). It is possible to link SPH nodes with standard finite elements such that solutions can be obtained for problems involving both highly distorted flow and structural response as mentioned by Johnson (1995) and . SPH method is well suited to impact problems such as bird-strike analyses because the particles are topologically independent from each other. It allows for severe distortions, eliminates many of the material interface problems, and being a Lagrangian technique SPH can be readily linked to standard FE formulations. SPH Algorithm for 2D-Axisymmetric Geometry A Lagrangian algorithm for SPH method is shown below as deduced from Johnson et al. (1996). The flowchart below shows the SPH algorithm, 10 American Institute of Aeronautics and Astronautics

11 Current velocities and Determine stresses from Determine stresses from strain displacements of all nodes. strains, strain rates etc. from and strain rate etc. SPH nodes Update velocities and Collect forces from standard Determine forces on each SPH displacements for all nodes. elements and SPH nodes node and its neighbors Figure 8 represents some features of the SPH technique. Node i is designated as the center node and the neighbor nodes are designated as nodes j. The distance between nodes is rij, the diameters of the nodes are di and dj, and the masses of the nodes are Mi and Mj. The masses remain constant throughout the computation, and are obtained from M = 0V0 where 0 and V0, represent the initial density of the material and the initial volume represented by the node. Figure 8: SPH Characteristics from Johnson et al. (1996), (with permission) The SPH approach allows for variable nodal connectivity, and this means that it is necessary for each center node i, to search and find the closest neighbor nodes j. This searching time can be significant, especially for larger problems. Smoothing Functions The smoothing function is an important part of the SPH algorithm. Two smoothing functions and their derivatives are shown in Figure 9. The B-Spline smoothing function is expressed as, (23) (24) Where and the smoothing distance is, (25) Where is dimensionless smoothing distance and di and dj are nodal diameters. Parameter is a user specified input in the range, , and its value is usually taken as 1. The nodal diameters are determined by, (26) where is the volumetric strain, d0 is the initial node diameter, and x0 and x are the initial and current x (radial) coordinates. A requirement of the smoothing function is that it exhibits the characteristics of a Dirac delta function 11 American Institute of Aeronautics and Astronautics

12 as hij approaches zero. The derivative of the smoothing function provides the weighing function for the strain rates and forces. The negative value of the B-spline derivative is, (27) (28) Figure 9: B-Spline and Quadratic smoothing functions and derivatives from Johnson et al. (1996), (with permission) The interesting feature of the B-Spline derivative is that exhibits a maximum at and decreases for both smaller and larger values of For it is logical that decreases because it has less influence as the distance is increased. To have decreased for is not intuitively satisfying because closer the neighbor node j comes to the center node i, the less influence it has. This can also lead to instabilities in compression. The quadratic smoothing function also shown in Fig 13 is expressed as, (29) The quadratic smoothing function derivative, the weighing function, always increases as the nodes move closer together, and always decreases as they move apart. This intuitively appears to be more realistic than the B-Spline derivative. It is also simpler than the B-Spline and does not have the compressive softening that can lead to compressive instabilities. Based on these arguments alone, it would appear that the newer quadratic smoothing function is an improvement over the traditional B-Spline. Strain Rates For axisymmetric geometry, the three normal strain rates, , the shear strain rate the rotational rate, and the volumetric strain rate, for the center node i, can be defined as, (30) (31) 12 American Institute of Aeronautics and Astronautics

13 (32) (33) (34) (35) Where, is the derivative of the smoothing function, Vj is the current volume at node j, are the x velocities of the node i and j, are the z velocities, and are the direction cosines from node i to node j, and xi is the x coordinate of node j. The volumetric strain is obtained by integrating the volumetric strain rate. The three factors are used to normalize the smoothing functions such that they will provide the exact strain rates in the three principal directions for states of constant strain rates. The resulting normalizing factors are, (36) (37) (38) The effect of using the Normalized Smoothing Function (NSF) algorithm can be significant. Fig.10 shows a cross- section of an axisymmetric ring of material represented by 25 SPH nodes (5 x 5). Figure 10: Equivalent strain rates in radially stretching cross-sections with a uniform arrangement of SPH nodes from Johnson et al. (1996), (with permission) The two cross-sections on the left side are for the Standard Smoothing Functions, which means that they do not use the NSF algorithm. The two on the right side do use the NSF algorithm. The top two use the B-Spline smoothing function and the bottom two use the Quadratic smoothing function. The numbers in the centers of the circular SPH 13 American Institute of Aeronautics and Astronautics

14 nodes represent the equivalent strain rate in the node when the axisymmetric cross-section is subjected to a stretching radial velocity of (39) The equivalent strain rate is defined as, (40) Substituting the resulting strain rates we get, (41) The upper left cross-section in Fig.10 (B-Spline without NSF) shows good accuracy in the center region ( = 1.01), but the side boundaries and corners are significantly low ( . The reason for the low strain rates at the boundaries is that the basic SPH algorithm is for an interior node, which assumes that there is a full distribution of neighbor nodes. When the distribution of neighbor nodes is not full, inaccuracies are introduced. The lower left in Fig.10 (Quadratic without NSF) has a similar pattern, but the interior nodes have a lower strain rate ( = 0.87). This does not necessarily mean that the B-Spline is more accurate than the Quadratic smoothing function, because the results in Fig. 10 are specifically for a uniform grid with a dimensionless smoothing distance of =1.0. For some other conditions (such as = 0.8) the Quadratic smoothing function is more accurate. At larger smoothing distances, the interior nodes approach an equivalent strain rate of =1.0, and at lower smoothing distances ( < 0.8) the strain rates decrease to at = 0.5.This effect is illustrated in Fig.11, where the equivalent strain rate of an interior node is shown as a function of the dimensionless smoothing distance, . Figure 11: Equivalent strain rate versus dimensionless smoothing distances for an interior SPH node in a radial stretching cross-section from Johnson et al. (1996), (with permission) Note that the B-Spline is more accurate at = 1.0 and the Quadratic is more accurate at = 0.8. The sharp slope changes for the Quadratic smoothing function at = 0.707, 1.000, 1.118, etc occur when the expanding smoothing distance encounters additional nodes at vij= 2.0.The Quadratic smoothing function derivative, , has a finite slope at vij= 2.0 (and therefore has a noticeable effect) whereas the slope of , for the B-Spline is zero at vij= 2.0.The NSF algorithm provides much improved results, with the ramp (at 0.5< < 0.6) caused by limiting the factors to 4.0.The NSF algorithm is also very effective when the nodes are in a non-uniform arrangement. Forces After the strain rates, rotational rate and volumetric strain rate are determined, it is possible to determine the shear and deviator stress, the pressure and the nodal artificial viscosity in the standard manner. These stresses must 14 American Institute of Aeronautics and Astronautics

15 then be converted to forces as shown in the flow chart previously mentioned. The nodal force in the x direction on node j due to the stress of node i, is, (42) The force due to the in plane stress is, (43) The force due to the hoop stress is, (44) And the force on the z direction is, (45) The forces on the center node i, due to the stresses in node i, are equal and opposite to the in plane forces in the x and z direction, (46) (47) Material Interface Fig 12 shows the material interface represented by SPH method. If the two materials are not bonded together, the standard SPH algorithm introduces unacceptably large errors because nodes from material A influence the strain rates in nodes in material B, and vice versa. Figure 12: Material interface for the SPH Method from Johnson et al. (1996), (with permission) Also, shear and tensile stresses are developed between the two materials such that sliding and separation are significantly inhibited. If SPH approaches are to be applied to problems involving sliding interfaces, then specialized interface algorithms must be developed. Linkage of SPH Nodes to Standard Finite Elements It has been noted previously that a desirable characteristic of the SPH approach is that it can be linked to standard finite element grids. This allows highly distorted materials to be considered together with structural response materials. Also, as shown in the algorithm flow chart, a single material model subroutine can be used for both the standard elements and the SPH nodes. This means that a material model needs to be incorporated and validated only one time, and that the user can be assured that the identical material model is used for both the standard elements and the SPH nodes. 15 American Institute of Aeronautics and Astronautics

16 Fig. 13 shows how SPH nodes can be attached to a standard finite element grid. The broken circle around node i represents the effective region ( ) for SPH nodes on node i .The size of the SPH interface nodes is equal to the size of the interface elements. The mass of the SPH interface nodes comes from the SPH nodes only, and this ensures that the mass on all nodes is essentially equal. Figure 13: SPH Node attachment to a standard Grid from Johnson et al. (1996), (with permission) The forces on SPH node i come from SPH nodes n 1, , n5 and from interface element B and C. There are no direct force contribution from standard nodes n 6, , n9, other than through the interface elements. Figure 14 shows the SPH nodes interacting with standard grid on a sliding interface. The SPH nodes are designated as slave nodes and the standard element from the master surface. As was the case for the attached interface, the standard nodes (m 2, m3 and m6) do not directly affect the forces on node i. The maximum allowable overlap of an SPH node and a slave node i , with master surface m2 - m3 , is shown as 0. When the SPH nodes are initially defined Figure 14: SPH nodes sliding on standard grid from Johnson et al. (1996), (with permission) as SPH nodes, then 0 = 0. When the SPH nodes are converted from the standard elements, 0 can be a small fraction of SPH node diameter. The approximations for the sliding interface are similar to those of the attached interface, but more research is required to evaluate and improve the existing algorithm. When slave node i overlaps ( ) master segment m2 - m3, then the three normal velocities of nodes i ,m2 and m3 are adjusted to : (1) conserve linear momentum, (2) conserve angular momentum and (3) provide a normal velocity match of node i on master segment m2 - m3 . The positions of the nodes are adjusted to be consistent with the velocity changes. 16 American Institute of Aeronautics and Astronautics

17 Acknowledgments The authors wish to acknowledge financial support from Pratt and Whitney. Also acknowledge support from Institute of Critical technology and Applied Science (ICTAS). Also thanking Wesley C.H Slemp for the assistance with the research. References 1 Sharing the skies manual. TP13549. Transport Canada; 2004. 2 Georgiadis S, Gunnion A J, Thomson R S and Cartwright B K, Bird Strike Simulation for Certification of the Boeing 787 Composite Movable Trailing edge, Compos Structure 86 (2008), pp. (258-267). 3 Abrate S. Impact on composite structures. Cambridge, UK: Cambridge University Press; 1998. 4 Roberts, G.D, Perieira M, Revilock, Binienda, W.K, Xie Ming and Braley, M, Ballistic Impact of braided Composites with a Soft Projectile, Journal of Aerospace Engineering, January, 2005, pp.3-7. 5 Talreja R, editor. Damage mechanics of composite materials, Composite materials series, vol. 9. Amsterdam: Elsevier; 1994. 6 Soden P.D, Kaddour A.S and Hinton M.J. Recommendations for designers and researchers resulting from the world-wide failure exercise. Compos Sci Technol 2004;64:589-604. 7 Puck A, Schurmann H. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 1998;58:1045-67. 8 Puck A, Schurmann H. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 2002;62:1633-62. 9 Hashin, Z. Failure criteria for unidirectional fiber composites. J Appl Mech 1980;47:329-34. 10 Johnson, A.F and Holzapfel, M, Modeling soft body impact on composite structures, Composite Structures 61 (2003), pp. 103113. 11 McCarthy M.A, Xiao J.R, Petrinic N, Kamoulakos A. and Melito,V, Modelling of bird-strike on an aircraft wing leading edge made from fibre metal laminates Part 1: material modeling, Appl Composite Mater 11 (2004), pp. 295315. 12 McCarthy M.A, Xiao J.R, Petrinic N, Kamoulakos A. and Melito,V, Modelling of bird-strike on an aircraft wing leading edge made from fibre metal laminates Part 2: modeling of impact with SPH Bird model, Appl Composite Mater 11 (2004), pp. 295315. 13 Meguid S.A, FE analysis of geometry effects of an articial bird striking an aeroengine fan blade, Int J Impact Eng (2007), doi:10.1016/j.ijimpeng.2007.04.008. 14 Hiermaier S, Thoma K, Computational simulation of high velocity impact situations using smoothed particle hydrodynamics, 9th DYMAT Technical Conference on Materials and Structural Modeling in Collision Research, TU Munich, Munich,Germany, 1995. 15 Wilbeck, J.S. Impact behaviour of lowstrength projectiles. AirForce Materials Laboratory, Technical Report AFML-TR-77 134, 1977. 16 Iannucci L,Bird strike impact modeling. I. Mech. E. Seminar 1998, Foreign Object Impact and Energy Absorbing Structures, Inst. of Mech. Engineers, London. 17 Roberts, G.D, Perieira M, Revilock, Binienda, W.K, Xie Ming and Braley, M, Ballistic Impact of braided Composites with a Soft Projectile, Journal of Aerospace Engineering, January, 2005, pp.3-7. 18 Kermanidis, T.H, Labeas, G, Sunaric, M and Ubles,L,Development and validation of a novel bird strike resistant composite leading edge, Applied Composite Materials, Vol. 12, No. 6. (November 2005), pp. 327-353. 19 Ubels, L. C., Johnson, A. F., Gallard, J. P. and Sunaric, M., Design and Testing of a Composite Bird Strike Resistant Leading Edge, SAMPE Europe: 24th International Conference and Forum, Paris, 13 April 2003, pp. 6. 20 PAM-SHOCK/PAM-CRASH FE Code. Engineering Systems International, F-94578 Rungis Cedex, France. 21 Nguyen, M,Q, A Review of Finite Element Software for Composite Impact Analysis, Journal of Composite Material, Vol 39 (2004), pp. 375-386. 22 Groenenboom P.H.L. Numerical simulation of hypervelocity impact using the SPH option in PAM-SHOCK. Int J Impact Eng 1997;20:30923. 23 ESI Group. Retrieved April 1, 2009. From http://www.esi-group.com/ 24 Johnson G.R, Stryk R.A and Beissel S.R, SPH for high velocity impact computations, Computational Method Appl Mech Eng 139 (1996), pp. 347373. 25 Monaghan J.J, Smoothed particle hydrodynamics, Ann Rev Astro Astrophys 1992;30:54374. 26 Lucy L.B. A numerical approach to the testing of fusion process, The Astron. J. 88 (1977) 1013-1024. 27 Gingold, R.A and Monaghan, J.J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Mthly Notices Roy. Astron. Sot.181 (1977) 375-389. 28 Benz, W. Smooth particle hydrodynamics: A review, Harvard-Smithsonian Center for Astrophysics, Preprint 2884. 1989. 29 Libersky, L.D and Petschek, A.D. Smooth particle hydrodynamics with strength of materials. Advances in the Free Lagrange Method, Lecture Notes in Physics 395 (1990) 248-257. 17 American Institute of Aeronautics and Astronautics

18 30 Johnson, G.R, Petersen, E.H and Stryk, R.A. Incorporation of an SPH option in the EPIC code for a wide range of high velocity impact computations, Int. J. Impact Engrg. 14 (1993) 385-394. 31 Petschek, A.G. and Libersky, L.D. Cylindrical smoothed particle hydrodynamics, J. Comput. Phys. 109 (1993) 76-83. 32 Johnson, G.R, Beissel, S.R. Normalized smoothing functions for SPH computations. International Journal of Numerical Methods in Engineering, 39, 2725-2741(1996). 33 Johnson, G.R. Linking of Lagrangian particle methods to standard finite element methods for high velocity impact computations, Nucl. Engrg. Des. 150 (1994) 265-274. 34 Attaway, S.W, Heinstein, M.W and Swegle, J.W. Coupling of smooth particle hydrodynamics with the finite element method, Nucl. Engrg. Des.150 (1994) 199-205. 18 American Institute of Aeronautics and Astronautics

Load More