Rubber

ARTICLE IN PRESS

International Journal of Impact Engineering 34 (2007) 647–667 www.elsevier.com/locate/ijimpeng

Impact of aircraft rubber tyre fragments on aluminium alloy plates: II—Numerical simulation using LS-DYNA D. Karagiozovaa, R.A.W. Minesb, a

Institute of Mechanics, Bulgarian Academy of Sciences, Acad. G. Bonchev Street, Block 4, Sofia 1113, Bulgaria b Impact Research Centre, University of Liverpool, Brownlow Street, Liverpool, L69 3GH, England Received 16 November 2005; received in revised form 13 February 2006; accepted 17 February 2006 Available online 8 June 2006

Abstract A discrete model for a reinforced rubber-like material is proposed in order to simulate numerically a debris tyre impact on a typical structure of an aircraft when using the FE code LS-DYNA. The model is calibrated using the static and dynamic test data for the actual tyre material. The dynamics of the tyre projectile is validated when comparing the numerical predictions with the response of a square aluminium alloy plate subject to a ribbon projectile impact having different initial velocities and impacting the plate at an angle of 301. Good agreement is obtained in terms of the strains in the plate caused by the ribbon impact. The numerically predicted deformations of the projectile also represent well the dynamics of the tyre ribbon recorded during the experiments. Some characteristic features of a soft projectile impact on a deformable plate are discussed. The developed model is then extended to the simulation of a full-scale impact test, and good agreement is shown between experiment and simulation. r 2006 Elsevier Ltd. All rights reserved. Keywords: Reinforced rubber; Discrete material model; Soft impact; Numerical simulation

1. Introduction The rebound of a projectile from a surface has been intensively studied in connection with the ballistic mechanics of ricochet for rigid and deformable bodies [1,2]. A new insight into the response of structural components associated with elastic–plastic large deformations is provided nowadays by the numerical simulations (e.g. [3–5]), aiming at an improvement of design and of an increase in product safety. The vast majority of studies on projectile rebound from a structure consider both the structure and the projectile to have a comparable stiffness and hence the subsequent recovery of the stored elastic energy causes motion changes in the projectile. In the case of an impact of a rubber like projectile on a metal plate, a significant proportion of the initial kinetic energy participates in the deformation of the projectile, so that the kinetic energy transferred to the plate depends significantly on the interaction between the plate and the Corresponding author. Tel.: +44 151 794 4819; fax: +44 151 794 4848.

E-mail address: [email protected] (R.A.W. Mines). 0734-743X/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijimpeng.2006.02.004

ARTICLE IN PRESS 648

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

Nomenclature A A1, A2 C ij D dy Ec E fz Ii Jel J K0 t Tk Tc W b dij li m0 nr, nc r0 rr, rc s

cross-sectional area of the fibres material parameters for the Mooney–Rivlin law (rubber part), Eq. (1) contravariant strain tensor material constants distance between reinforcing cords in the y direction elastic modulus of the fibre material elastic modulus of the aluminium alloy number of cord layers per unit thickness strain invariants elastic volume ratio total volume ratio bulk modulus true stress kinetic energy (plate or tyre projectile) deformation energy (plate or tyre projectile) strain energy potential scaling coefficient, Eq. (6) covariant strain tensor deviatoric extension ratios shear modulus the Poisson ratios for the rubber part and reinforcing fibres, respectively actual density of the tyre material density of the rubber part and reinforcing fibres, respectively nominal stress

projectile, and on the dynamics of the projectile during impact. Different angles of impact result in different projectile trajectories, which also contribute to the variation of the forces and energy transmitted to the plate. The aircraft tyre materials used in practice have a complex structure consisting of a rubber-like part and reinforcement [6]. The reinforcing cords are ropes twisted from Nylon fibre and their stiffness in compression is significantly lower then the stiffness in tension, moreover they are placed at different angles [6,7]. This material structure causes anisotropic properties in the different directions of loading and non-linear stressstrain dependence when subjected to large deformations. The complex material structure requires an adequate modelling, which can represent correctly the stiffness of the impactor since this characteristic can significantly influence the interaction with a structure due to the large flexibility of the soft projectile. Several approaches are usually used in the finite element analysis of cord-reinforced composites. One class of computational models are based on the analysis techniques for laminated fibrous composites [8,9]. The material properties of the individual components are averaged over a layer of a laminated shell element and described by an anisotropic material law. This approach however, requires assumptions for a number of material constants related to the ‘averaged’ layers, which causes difficulties prior to the analysis and after, for the recovering of the properties of the individual components. This approach could have restricted application in the case of very large anisotropy, which might result in a non-positive defined elasticity matrix. An alternative approach is to use one-dimensional bar elements to simulate the behaviour of the reinforcing cords. These elements are superimposed on the rubber elements by satisfying the compatibility conditions at their common nodes. This approach was first proposed for reinforced concrete [10] and applied later to rubber composites [7,11]. The concept is simple and allows the modelling of anisotropic material properties of reinforced rubber when varying the bar position and orientation of the bar elements. The major disadvantage of this approach is that the locations and orientation of the cords are restricted by the finite element discretisation, which could present difficulties for a high number of reinforcing cords. However, the low

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

649

computational cost and relatively small number of the required material constants make this approach attractive, particularly when the material parameters need to be obtained by calibration. Another efficient approach has been developed to model cord-reinforced composites, which is characterised by superimposing so-called ‘rebar’ elements, consisting of one or more reinforced cord layers with arbitrary orientations, on corresponding rubber elements. This technique was originally applied for reinforced concrete [12] and is widely used nowadays to model cord reinforced rubber [13–16]. Two and three dimensional rebar elements have been implemented in FE codes such as MARC and ABAQUS. This approach is discretisation independent and is probably the most effective for modelling reinforced composites. A further step in modelling reinforced rubber-like materials is related to the development of a mixed twofield displacement-pressure energy function, which leads to constitutive equations applicable for fibre reinforced materials that experience finite strains [17,18]. Obviously, this approach requires an implementation of new material model in the FE codes. The primary objective of the present study is to develop a reliable and inexpensive model for a tyre material when using the available material library in LS-DYNA and to estimate the applicability of such model for simulating an impact on a deformable structure. Given the complexities of the problem, viz., large dynamic strains and contact conditions between the missile and the plate, it was decided to take the simplest possible approach within the limitations of DYNA. Hence the second approach was used, i.e., one-dimensional bar elements representing the reinforcing cords. The proposed discrete model for the tyre material is calibrated using the static and dynamic test data for the actual tyre material [6]. The validation of the model is performed when comparing the numerical predictions with the response of a square aluminium alloy plate to a tyre projectile impact [6] having different initial velocities and impacting the plate at 301. 2. Material characteristics The tested tyre material has a complex structure consisting of a rubber-like material reinforced with nylon cords. The cords are placed at different angles, but in general these angles are smaller than 301 with respect to the circumferential axis of the tyre (Fig. 1(a)). The reinforcing cords are ropes twisted from Nylon fibres and their stiffness in compression is significantly lower then the stiffness in tension. This material structure results in anisotropic properties and a non-linear behaviour when subjected to large deformations, as revealed by the test data presented in [6]. The experimental static material properties in tension in the x and y-directions are shown in Fig. 1(b) and (c), respectively. The x-direction for the tested and modelled specimens is associated with the circumferential direction of an actual tyre and the y-direction corresponds to the tyre radial direction (Fig. 1(a)). The z-direction is the through thickness direction for an actual tyre. During a debris impact on a structure, large compressive deformations of the projectile and large bending deformations, and therefore considerable tensile strains can occur depending on the angle of impact due to the flexibility of the projectile. Thus, to accurately model the way a tyre projectile deforms upon impact, both tension and compression data are required. 3. Discrete model for a tyre material The approach, which uses one-dimensional bar elements to simulate the behaviour of the reinforcing cords in rubber, is used to model the tyre material. In general, the properties of the rubber-like part determine the flexibility of the material, whilst reinforcing cords are used to control the anisotropic properties in tension in the x and y-direction. It is assumed that there is no anisotropy in the z-direction but that the model material retains the orthotropic properties in tension in x and y-directions, similar to the actual material. These material properties are modelled using layers of cords placed at (x, y) planes, which are situated at a constant distance along the z-axis as shown in Fig. 2(a). The modelled specimens have dimensions 75  20  10 mm3 and 20  75  10 mm3 for tension in the x and y-direction, respectively, which are the same as the specimen dimensions from [6]. The cords are placed at a certain constant angle with respect to the x-axis, which is determined by the uniform finite element discretisation for the rubber part of the tyre material. It is assumed that the cords are

ARTICLE IN PRESS 650

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

Z y

x

Tread form

(a)

Fibre laminates

Sp C

Stress (MPa)

80 60

Sp B 40 Model 20 Sp D 0 0

20

(b)

60 40 80 Strain (75mm gauge length) (%)

100

Stress (MPa)

80 60 40

Sp E Sp G

20 Model

Sp F

0 0

(c)

20

40 60 80 Strain (75mm gauge length) (%)

100

Fig. 1. Comparison in terms of nominal stress-strain between the experimental results and numerical predictions; (a) orientations with respect to tyre, (b) tension in the x-direction (Model B), (c) tension in the y-direction (Model A).

incompressible but have linear elastic properties at tension, so that a Discrete beam formulation with Cable material properties was used [19]. It should be noted that later versions of DYNA allow for the non linear modelling of the cords. These cable elements are superimposed on the solid elements modelling the rubber part by satisfying the compatibility conditions at their common nodes, using the definition associated with the LS-DYNA concept Constrained_Lagrange_in_solid [19]. The rubber part of the tyre is modelled using a Mooney–Rivlin material model [19–21] together with fully integrated quadratic 8-nodes solid elements with nodal rotations. It was established that this type of solid element has better performance when large strains at the common nodes between the solid and discrete elements occur and when large gradient contact forces develop during the plate impact. No strain-rate effects are assumed for the material model although an increase of the magnitudes of the resultant dynamic forces at compression has been observed [6]. This increase in forces is associated with inertia effects in the tyre material due to the size of the specimen rather than with strain-rate sensitivity. Due to the overlapping of the rubber material and the cords, a reduced material density for the rubber-like part and cords is used, in order to obtain the actual density of the tyre material r0 ¼ 1020 kg m3 . It is assumed that both material components have equal density, which results in a reduced model density r ¼ r0 =ð1 þ rÞ, where r ¼ V fibres =V rubber is the reinforcement fraction represented by the ratio between the volumes of the rubber and cords components. Due to the high ratio of the bulk modulus to the shear modulus, a refined mesh is recommended [22] in order to overcome the numerical difficulties resulting from the incompressibility of the rubber, so that an

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

651

Fig. 2. Models for the tensile specimens; (a) geometry of the fibre layers, (b) tension in the x-direction, (b) tension in the y-direction. Cords are given reduced diameter for clarity.

Fig. 3. Models for the specimen compressed in the z-direction; (a) an actual cylindrical specimen, (b) an equivalent reinforced square prism—side view, (c) square prism—top view.

element size of 2.22  2.27  1.67 mm3 was used when modelling the tensile material properties. Another reason to use small size elements in the z-direction was imposed by the possibility to model materials with high reinforcement fraction, which is typical for the tyre materials. The particular element sizes throughout the analysis ware selected after a careful mesh sensitivity analysis on the performance of specimens’ models made of a rubber-like material with and without reinforcement. In the former case, the mesh sensitivity analysis was performed for a constant angle of reinforcement. The present approach applied to model the reinforced material, although inexpensive computationally, created problems when modelling the cylinder specimen used in compression tests [6]. In order to overcome the modelling problem and also to estimate the influence of the reinforcement on the compressive force, the following procedure is applied. The actual cylindrical specimen having a diameter 25 mm [6] pffiffiffi is replaced by an equivalent prism with a square base having c ¼ R p ¼ 22:15 mm and height 18 mm as shown in Fig. 3. The selected solid element size is 2.46  2.46  1.67 mm3, which determines the reinforcement angle as 18.41. The proposed discrete model for the actual tyre material is characterised by a number of constants, namely the Poisson ratio, nr, the material constants approximating the rubber-like material and the cords characteristics: the cross-sectional area A, Young modulus, E, Poisson ratio nc and angle of reinforcement.

ARTICLE IN PRESS 652

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

4. Calibration with the static tests The above mentioned constants, which characterise the rubber and reinforcing part of the model, are obtained by calibration using the test results for the tyre material from [6]. As discussed in [6], due to the lack of availability of separate data for the Nylon reinforcement and the rubber matrix, these properties were derived from testing the reinforced tyre rubber in various directions. The rubber matrix of the tyre was modelled using the Mooney–Rivlin model in DYNA as described in [6]. The theory given in [6] was brief, and the more general three-dimensional formulation is given here, in order to fully understand the numerical simulations and the limitations of the Mooney–Rivlin model. The Mooney–Rivlin model is derived using finite strain elasticity, which is discussed in detail by Treloar [20] and Green and Adkins [21]. In general terms, the strain energy potential, W, is given by 1 el ðJ  1Þ2 , (1) D1 where A1, A2, and D1 are material parameters, I¯ 1 and I¯ 2 are the first and second deviatoric strain invariants defined by W ¼ A1 ðI 1  3Þ þ A2 ðI 2  3Þ þ

2 2 2 I¯ 1 ¼ l¯ 1 þ l¯ 2 þ l¯ 3

and

ð2Þ ð2Þ ð2Þ I¯ 2 ¼ l¯ 1 þ l¯ 2 þ l¯ 3 ,

(2a,b)

where the deviatoric extension ratios are given by l¯ i ¼ J 1=3 li

(2c) el

J is the total volume ratio, J is the elastic volume ratio and li are principal stretches. The initial shear modulus and bulk modulus are given by m0 ¼ 2ðA1 þ A2 Þ and

K 0 ¼ 2=D1 .

For the particular case of simple extension, the true stress t is given by [19–21]       1 qW 1 qW 1 1 2 2 t ¼ 2 l1  þ ¼ 2 l1  A1 þ A2 l1 qI 1 l qI 2 l1 l

(3a,b)

(4)

which is consistent with Eq. (4) in [6]. Once constants A1 and A2 are obtained from uniaxial tests, threedimensional true stress is derived using [19–21]:   pffiffiffiffiffi qW 2 qW 2 qW 2 qW tij ¼ 2 I 3 dij þ pffiffiffiffiffi þ I 1 pffiffiffiffiffi C ik C jj , (5) C ij  pffiffiffiffiffi qI 3 qI qI I3 1 I2 1 I 2 qI 1 where d and Cij are strain tensors. It can therefore be seen that the derivation of three-dimensional stresses and strains from the model that has been calibrated using two-dimensional stress system, is complex, and includes many assumptions. Care is required to differentiate between nominal stress and true stress. For more information see [20,21]. Specifically, the Mooney–Rivlin model can become inaccurate for strains over 70%, and for cases where shear and biaxial effects dominate. It should be noted in passing that Green and Atkins [21] consider the case of inextensible cords in a finite element continuum (i.e., reinforced tyre rubber). They provide analytic solutions for flexure of reinforced cuboids and for in plane loading of reinforced membranes. Such analytic solutions could be used to further validate and investigate the numerical models used here. The static tensile material tests were simulated when a constant low velocity was applied to both ends of the specimens and the displacements in the direction perpendicular to tension were constrained. It is evident in Fig. 4 that the tensile properties of the reinforcing fibres (the reinforcing angle is 181) contribute significantly to the response of the model to tension in the x-direction. Due to the small angle with respect to the x-axis, the cords could be subjected to considerable tensile deformations. In contrast, the reinforcing cords contribute negligibly to the tensile force in the y-direction as they are subjected mainly to rotation until considerably large tensile strains occur. In fact, the reinforcing cords do not contribute at all to the tensile force in the y-direction up to 110% strain, which is the maximum tensile strain from the test result in the y-direction (Fig. 5).

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

653

Fig. 4. Shapes of the tensile specimen in the x-direction from simulation (Model B); (a) initial, (b) deformed at 35% strain. Cords are given reduced diameter for clarity.

Fig. 5. Shapes of the tensile specimen in the y-direction from simulation (Model B); (a) initial, (b) deformed at 110% strain. Cords are given reduced diameter for clarity.

A comparison between the stresses developed in the rubber part of the tyre material (recorded in the middle cross-section of the modelled specimen) and the axial stress that occurs in a single reinforcing cord is shown in Fig. 6(a). One can see that the tensile material strength in the x-direction is mainly due to the reinforcement. By way of contrast, the tensile strength in the y-direction is determined only by the tensile properties of the rubber material (at least for strains smaller than 110%). It is shown in Fig. 6(b) that the tensile stresses in the y-direction are equal to the stresses in the rubber part of a specimen subjected to tension in the x-direction. Fig. 7 shows compression stress strain experimental data (using a cylindrical specimen) [6] for two cases, namely (a) core only and (b) through thickness material. It was observed in the experiments that the compressive characteristics of the core-only material (i.e. reinforced region) and the through the thickness (i.e. reinforced region and rubber outer tread) response differ considerably within the expected range of the compressive strains and the core only material manifests larger stresses. There is a range of compressive strains between 50 and 75%, approximately for the core-only material, where the stresses vary only slightly with the strains, which can be attributed to the structural response of the specimens (e.g. local damage) but not to the material characteristics themselves. It is anticipated that this response is due to the physical interaction between the fibres in the actual material, which is difficult to model using the available elements and material models in LS-DYNA. For that reason, in terms of compressive strains, it is aimed to obtain material constants for the model, which can predict nominal stress strain relationship closer to the one characterising the through-the-thickness material for strains smaller than 70–75%, as this material is closest to the characteristics of the actual projectile. The nominal stress in the simulation is obtained when the reaction force at the edge of the tested specimen is divided by its initial area. Consequently, the stress-strain relationship in tension in the y-direction (Fig. 1(c)) and the ‘through the thickness’ compression stress-strain curve (Fig. 7(b)) are used to obtain the material constants for the rubberlike material when the Poisson’s ratio is assumed as nr ¼ 0:495. The two parameters for the Mooney-Rivlin material model can be obtained, e.g., either graphically or using the curve-fitting procedure implemented in

ARTICLE IN PRESS 654

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667 100

Stress (MPa)

80

Beam (axial stress) 60

40

20

Rubber part (von Mises stress)

0 0

5

10

15

(a)

20

25

30

35

40

Strain (%)

Von Mises stress (MPa)

25

20

15

y-direction 10

x-direction

5

0 0

(b)

20

40

60

80

100

Strain (%)

Fig. 6. Stresses resulting from quasi-static tension in the x and y-direction; (a) comparison between stresses in the rubber part and a cable (tension in the x-direction), (b) comparison between stresses (using Model A) in the rubber parts of the tensile specimens in the x and ydirection.

LS-DYNA. This latter procedure led to A1 o0 and A2 40 for the Mooney-Rivlin law, when using the standard curve-fitting procedure in LS-DYNA. In order to satisfy the Drucker postulate (ds=dl40, where s and l are the uniaxial stress and stretch, respectively), the negative constant should be assumed equal to zero. However, the obtained stress-strain curve using only one constant resulted in a good agreement with the experimental curve only for small compressive strains and there was large deviation for both large compressive and large tensile strains. Nevertheless, the constants A1 ¼ 0 and A2 ¼ 6:8 N mm2 obtained by the standard fitting procedure for the Mooney-Rivlin model are used as an initial approximation for the material model with no cords (Model A). This is required as a starting point for curve fitting, when the cords are introduced into the model. The nominal stresses obtained from model A are compared with the tensile test data in the y-direction (Fig. 1(c)) and the static compression stress (Fig. 7(b)). It should be noted that the simulation results in Fig. 7 are obtained using an equivalent cuboid specimen model as discussed previously and shown in Fig. 3. Next, the cords are introduced (Model B). The constants calibrating the cords’ characteristics were obtained when ‘fitting’ the tensile curve to the stress-strain relationship from the tensile test in the x-direction. In this case, the Poisson ratio for the cord material is assumed as nc ¼ 0:3 but the cross-sectional area and elastic modulus of the cords are varied. The reinforcement in the models for the tensile tests consists of four layers of cables (see Fig. 2). The large anisotropy in tension is achieved by increasing the axial forces in the cables and varying the angle of reinforcement. An elastic modulus of the cords of 250 MPa was assumed [6], and this gave an individual cross-sectional area A ¼ 3 mm2 which results in a cord diameter of 1.95 mm. Due to the linear elastic properties of the cables and the major contribution of this component to the tensile characteristics in the x-direction, only an almost linear approximation of the tensile stress-strain relationship in this direction is possible as evident in Fig. 1(b). The resulting volume fraction of the cords’ component is quite high as V fibre =V rubber ¼ 0:567, which determines a reduced material density as r ¼ 0:67r0 , where r0 is the density of the actual tyre material.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

655

500 BC1

Stress (MPa)

400 300

BC3 200 100 Model 0 0

20

(a)

40 60 Strain (%)

80

100

500

Stress (MPa)

400 BC3

300

BC1 200 Model 100 0 0

(b)

20

40 60 Strain (%)

80

100

Fig. 7. Comparison between the experimental results and numerical predictions (Model A) for compression in the z-direction; (a) coreonly material, (b) through-the-thickness material.

200

Stress (MPa)

160

Cylinder 120

Prism

80 40 0 0

0.1

0.2

0.3

0.4

Strain

0.5

0.6

0.7

0.8

Fig. 8. Comparison between the static compression nominal stress-strain curves for a cylinder and an equivalent square prism (Model A).

A comparison between the static compression of a model for a non-reinforced actual cylindrical specimen and a rubber only equivalent prism is made in order to estimate the possible error when replacing the actual cylindrical specimen by a square prism (the compression of a square prism is not entirely uniaxial as it is supposed to be in the experimental set up) Fig. 8. A square prism could be used to represent the actual cylinder for strains smaller than 70%. It should be also noted that the reinforcement of the model does not contribute considerably to the compression in the z-direction and only slightly increases the compressive stress due to the variation of the Poisson effect. The material constants for the model are obtained by minimising the overall deviation of the numerical predictions from the material test results [6]. The smallest deviation of the predicted nominal stresses from the experimental tensile stresses in the x-direction (Fig. 1(b)) and compressive one (Fig. 7(b)) are obtained when

ARTICLE IN PRESS 656

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

using A1 ¼ 0.5 N mm2 and A2 ¼ 4.2 N mm2 for the Mooney–Rivlin material model and reinforcement at 181 with respect to the x-axis with cords having diameter of 1.95 mm and stiffness of 250 MPa. It should be noted that the numerical simulations for the calibration of the static stress-strain relationship are done at a constant velocity of 0.5 m s1 in compression and tension due to the very significant increase of the computational time when decreasing the load speed. Only the tensile test in the x-direction is simulated with a constant velocity of 0.1 m s1 as well, but virtually no difference is observed in the response in comparison to the tension case with 0.5 m s1.

5. Verification with the dynamic tests Although the tested material has not shown significant inertia effects during the dynamic tests [6], it is important to verify the proposed model with the dynamic test results due to the large flexibility of the tyre projectile. The dynamic verification of the model for the tyre material is done simulating impacts in the z-direction when using the equivalent square prism model shown in Fig. 3. The reduced density for the rubber part of the reinforced prism is r ¼ 0:737r0 . The loads are applied as impacts by a mass G ¼ 12 kg with initial velocity 7 m s1 and by a mass G ¼ 0.385 kg with initial velocity 50 m s1. No special surface treatment was applied to the contact surfaces of the tested specimens, so that the measured static friction coefficient of 0.778 [6] is used for the simulations. Comparisons between the numerical predictions and the test results are presented in Fig. 9(a) and (b) for a low and high impact velocity, respectively. It is evident that the selected specimen geometry, the model for the reinforcement and the material characteristics of the tyre components lead to numerical predictions for the dynamic load–displacement histories, which are in good agreement with the corresponding test results. A discussion of the strain rates for these cases is given in [6].

100 Model

Load (kN)

80 60 40 Drop Hammer BT2

Drop Hammer BT1

20 0 0

2

4

(a)

6 8 Displacement (mm)

10

12

14

120 Model

Load (kN)

100 80 Gas Gun BT2 60 40

Gas Gun BT3 20 0 0

(b)

2

4

6 8 10 Displacement (mm)

12

14

Fig. 9. Comparison between the experimental results and numerical predictions (Model B) for the dynamic compression of a through-thethickness specimen (a) impact with V 0 ¼ 7 m s1, G ¼ 12 kg, (b) impact with V 0 ¼ 50 m s1, G ¼ 0:385 kg.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

657

140 120

Dynamic 50m/s

Load (kN)

100 80 Dynamic 7m/s

60 40

Static

20 0 0

2

4

6 8 10 Displacement (mm)

12

14

16

Fig. 10. Summary of the static and dynamic force-displacement histories for Model B and compression in the z-direction of a through-thethickness specimen.

Table 1 Characteristics of the tyre material model Rubber part Mooney–Rivlin law

Reinforcement Cable material

A1 (N/mm2)

A2 (N/mm2)

nr

rr (kg/m3)

E (MPa)

A (mm2)

nc

rc (kg/m3)

Angle

0.5

4.2

0.495

0.7r0–0.75r0

250

3

0.3

0.7r0–0.75r0

E181

Although no strain rate effects are taken into account, an increase of the dynamic force occurs when simulating the dynamic tests due to the inertia effects in both rubber part and the reinforcing cables, which are related to the size of the material sample selected for simulation. A summary of the simulated static and the two dynamic load-displacement histories (z-direction) in Fig. 10 gives an estimate for the dynamic effect in the modelled specimens. It should be noted that a similar increase of the compressive load is observed experimentally when increasing the impact velocity [6]. 6. Summary of the model characteristics and model parameters The various comparisons for the static stress-strain relationships (Figs. 1(b), (c), 7(a) and (b)) and dynamic force–displacements histories (Fig. 9) confirm that the selected approach for modelling a reinforced rubberlike material together with the calibrated material constants is capable to reproduce the response of the actual tyre material to standard tensile and compressive tests. The constants for the material model are summarised in Table 1. 7. Validation of the material model simulating a tyre impact on an aluminium alloy plate 7.1. Experimental arrangement [6] In order to verify the dynamics of the tyre projectile and estimate the reliability of the selected model for the tyre material, an impact of a tyre projectile on an aluminium alloy plate is simulated. The experimental method, data reduction and results are given in [6]. The numerically predicted strains in the plate, which result from various impacts, are compared with the measured ones in four particular positions at the distal surface of the plate and particular results are presented here for impacts at 301. The test arrangement for the tyre debris impact on a square aluminium alloy plate is described in [6]. An aluminium alloy plate 300  300  1.6 mm3 was bolted on a thick steel frame and positioned at different angles with respect to the axes of the projectile. The inside dimension of the frame is 260  260 mm2, which led to an

ARTICLE IN PRESS 658

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

actual deformable size of the plate 260  260  1.6 mm3. The projectile is cut from a tyre in a shape of a rectangular prism (ribbon) 60 mm long  15 mm wide  30 mm deep having a weight of 28 g. In all experiments, the plate centre (x ¼ 0, y ¼ 0) coincides with the centre of the ribbon (x ¼ 0, y ¼ 0) and a ribbon projectile is accelerated up to 136 m s1 at an (0, 0, 301) orientation in rectangular co-ordinate system (x, y, z). Strain gauges are attached at 01 and 901 to the distal side of the plate in the area of impact and at a distance 100 mm from the edge of the plate in the direction of impact. The records from these strain gauges [6] are used to verify the numerically obtained strains in the x and y-directions. 7.2. Finite element model 7.2.1. Geometry, boundary and contact conditions The structure and the material properties of the ribbon are selected according to the material model calibration presented in Sections 2–6 and Table 1. Eight reinforcing layers are placed parallel to the plane (x, y) at equal distances along the z axis when using 18 fully integrated solid elements for the rubber part in this direction (Fig. 11). The reinforcement is placed at 18.41 with respect to the x-axis, which is related to the uniform discretisation (2.5  2.5 mm2) in the (x, y) plane. The elements are 1.65 mm in the z-direction. It should be noted that the latter dimension is similar to the diameter of the cords (1.95 mm). Belytschko-Tsai shell elements 6.5  6.5 mm2 are used to model the plate where clamped boundary conditions are assumed along all edges. The attempt to use larger shell elements with adaptive mesh increased significantly the computational time and caused difficulties in the contact area for the 301 impact. An important issue in the simulation of the rubber tyre fragment impacting an aluminium plate is the nature of the contact conditions between the two. This is complex given the large deformation in the rubber missile during impact and the dynamic nature of the contact. However, given the complexity of the problem, it was decided initially to use one of the simplest contact algorithms in DYNA, viz. Nodes_to_surface_contact [19]. It should be noted that computational contact mechanics is a complex subject [23] and that a number of approaches are possible. The theory used here relates to finite deformations within the contact surface, and to a node of the contact body (slave body) sliding over several elements of the other body (master surface). The model assumes Coulomb friction and hence the only parameter required is the coefficient of friction (m ¼ 0:778 from [6]). The high strain gradient in the model projectile during the initial contact imposed a condition on the time step and it was reduced by a factor of 0.5 compared to the default time step in LS-DYNA. 7.2.2. Material properties of plate and associated strains A piecewise linear stress-strain relationship in the true stress–true strain space is assumed for the aluminium alloy plate according to Table 2 while the elastic modulus is E ¼ 72 GPa [24]. No strain rate or plate failure effects are taken into account in the numerical simulations. The reasons for not modelling plate rupture were that large-scale plate rupture did not occur in experimental tests [6] and also the focus for this study was the development of a tyre rubber impact model and not a detailed model of plate behaviour. The calculated strains in the plate are associated with the two elements on the right-hand side of the plate centre and the two elements located at 100.75 mm from the plate edge in the right–hand-side half of the plate

Fig. 11. Geometry of the ribbon projectile; (a) geometry of the reinforcing layers, (b) top view. Cords are given reduced diameter for clarity.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

659

Table 2 True stress-true strain relationship for the aluminium alloy [24] Effective stress (N/mm2) Plastic strain

276.0 0.0

315.4 0.002405

332.8 0.00426

351.1 0.00752

380.7 0.01781

462.2 0.06635

495.0 0.08906

Fig. 12. Location of the elements in the plate for the strain comparison (distal surface of the plate).

model. The locations of the elements, which are used for a comparison with the records from the strain gauges, are shown in Fig. 12. The strains for comparison in the upper elements are associated with the y-strains and correspond to the strain gauges SG2 (centre) and SG4 (at a distance). The elements below the x-axis are used to compare the x-strains being associated with the strain gauges SG1 (centre) and SG3 (at a distance).

7.3. Numerical simulation of the experiment The impacts at 301 by a ribbon projectile at an initial velocity of 135 m s1 are characterised by large bending deformations of the projectile (Fig. 13(a)) and considerable sliding along the plate, when it leaves slide marks on the plate surface at the location of contact (Fig. 13(b)). For further information see [6]. The response of the plate and the ribbon is shown in Fig. 14 for an impact velocity of 135 m s1 at several particular times. It is evident, that the initial deformation phase of the projectile is characterised by bending and sliding along the plate and this response lasts until t ¼ 0:6 ms. Due to the flexibility of the ribbon, a second contact between the ribbon and plate occurs at t ¼ 0:68 ms and again the projectile slides along the plate a certain distance before rebounding. Loss of contact occurs at 1.16 ms. The two sliding contacts observed from the numerical simulations have their experimental evidence as slide marks on the plate as shown in Fig. 13(b) – one being associated with the area of the initial contact and another closer to the edge of the plate. Fig. 15 gives a three dimensional view at 0.68 and 0.8 ms, highlighting sliding behaviour. It can therefore be concluded that the numerical simulation models the experimental behaviour from the point of view of rubber missile deformation. Figs. 14 and 15 show that the tyre projectile undergoes large deformations and this behaviour results in an energy partition, which is different from the impact energy, associated with a hard projectile impact. In the present case of a ‘soft’ impact, a significant proportion of the initial kinetic energy is transformed into deformation energy of the projectile (70 J) and a relatively large energy is dissipated through the sliding contact (30 J). The variation of the kinetic energy of the plate (10 J) and the projectile and the variation of corresponding deformation energy are compared with the sliding energy in Fig. 16(a) and (b). One can see that

ARTICLE IN PRESS 660

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

Fig. 13. Ribbon impact at 301 with V 0 ¼ 135 m s1: (a) video from the impact, (b) rubber slide marks on the plate from [6].

Fig. 14. Response of an aluminium plate to a 301 impact V 0 ¼ 135 m s1 using Model B.

the sliding energy is considerably larger than the kinetic energy of the plate and is comparable to the plate deformation energy (20 J). The largest proportion of the initial kinetic energy (150 J) is transformed into deformation energy of the projectile. After loss of contact, the missile continues to vibrate and energy is

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

661

Fig. 15. Response of an aluminium plate to a 301 impact with an initial velocity 135 m s1 – secondary impact using Model B.

280 240

Energy (J)

200 160

Kinetic energy (Projectile)

120 80 40

Kinetic energy (Plate)

0 0

0.0005

0.001

(a)

0.0015

0.002

0.0025

0.003

0.0035

Time (sec) 120

Deformation energy (Projectile)

Energy (J)

90

60

Sliding energy 30

Deformation energy (Plate) 0 0

(b)

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

Time (sec)

Fig. 16. Partition of the initial kinetic energy for a 301 impact with an initial velocity 135 m s1: (a) kinetic energy, (b) deformation and sliding energy using Model B.

dissipated in material damping. It should be noted that the relative proportion of these energies will be dependent on missile geometry and impact conditions, and that the missile flexibility will ameliorate the energy input into the plate. A comparison between the experimentally recorded strain-time histories from the strain gauges SG1–SG4 [6] and the corresponding numerical predictions for a ribbon impact at 135 m s are shown in Fig. 17(a) and (b), respectively. The numerically predicted strain magnitudes are somewhat larger than the experimentally recorded ones, but this behaviour could be attributed to the highly non-homogeneous strain state at the initial impact area. However, the overall shapes of the curves are similar.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

662

0.016 0.014 0.012 0.01 SG2

Strain

0.008 0.006 0.004

SG1

0.002 0 SG3

-0.002 -0.004 -1

-0.5

0

0.5

1

1.5

2

2

2.5

3

Time (msecs)

(a) 0.016 0.014

SG2

0.012 0.01

Strain

0.008

SG1

0.006 0.004 0.002 SG3

0 -0.002

SG4

-0.004 0

(b)

0.5

1

1.5 Time (msecs)

Fig. 17. Variation of the strains in four particular locations in the plate for a 301 impact with V 0 ¼ 135 m s1: (a) experimentally obtained results (see [6]), (b) predictions from the numerical simulation using Model B.

The residual plastic strains from the simulation, which occur at the near and distal surface of the plate, are shown in Fig. 18. The largest plastic strains occur at the distal surface in the areas of the initial and secondary impacts. The locations of these plastic strains correspond to the slide marks observed experimentally (Fig. 13(b)). Plastic strains with smaller magnitudes occur at the near surface along the plate edge, which is closer to the secondary impact. 7.4. Interaction between the plate and projectile depending on the loading conditions The response of the plate to a 301 impacts with initial velocity of 135 m s1 suggests that the impact velocity and the friction coefficient are the two parameters which influence the dynamics of the projectile and therefore the energy partition associated with the plate and ribbon deformation as well. The influence of the impact

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

663

Fig. 18. Plastic strains in a plate resulting from a 301 impact with V 0 ¼ 135 m s1: (a) near surface, (b) distal surface using Model B.

velocity and the friction coefficient on the dynamics of the tyre projectile and the residual plastic strains in the plate are briefly discussed in this section. 7.4.1. Influence of the impact velocity The responses of the examined aluminium plate to 301 impacts with initial velocities 95, 75 and 55 m s1 were simulated and compared. The projectile tends to tumble along the plate rather than to slide when decreasing the impact velocity but nevertheless it hits the plate twice and plastic deformations occur in both areas of the initial and secondary impact. The location of the maximum plastic strains also changes when decreasing the impact velocity. The largest plastic strains occur in the area of the initial impact for initial velocities 95 m s1 and higher, while lower velocity impacts cause larger plastic strains in the area of the secondary impact. 7.4.2. Influence of the friction coefficient The response of the examined aluminium plate to a 301 impact with initial velocity 135 m s1 was also simulated and compared when assuming friction coefficients of 0.6, 0.778 (measured experimentally) and 0.95. The results from the numerical simulations show that the variation of the friction coefficient does not influence the dynamics of the projectile for the particular impact velocity. The increase of the friction coefficient however causes variation of the plastic strains in the plate without affecting their locations. The largest plastic strains (max(ep) ¼ 0.048) develop at the distal surface in the area of the initial impact for friction coefficient 0.778, while the plastic strains in the area of the secondary impact are only 0.02. Comparable plastic strains occur in the areas of the initial and secondary impacts when assuming friction coefficients of 0.6 or 0.95. Maximum plastic strains 0.028 and 0.023 are associated with the friction coefficient 0.6, while values of 0.033

ARTICLE IN PRESS 664

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

and 0.028 for the maximum plastic strains occur for the friction coefficient 0.95, but spread over a larger area. Plastic strains develop at the near surface in the area of the secondary impact as well, when assuming friction coefficient 0.95, which is not observed for friction coefficients 0.6 and 0.778. Presumably, the variation of the plastic strains is related to the different partition of the initial impact energy due to the large proportion of the total energy absorbed during the sliding contact. 8. Extension of proposed tyre model to larger tyre fragments The validation of the material model done by simulating a tyre debris impact on an aluminium alloy plate shows that the proposed discrete model is realistic and adequate to accommodate the anisotropy and nonlinearity inherent in the tyre problem. However, although inexpensive numerically, the proposed approach possesses two major disadvantages. In the real cord-rubber lamina the bonding between the surface of the cords and surface of the rubber is more than 60% [10], whereas in the present model, the cords are attached only in the common nodes. In addition, the diameter of the cords (1.95 mm) is similar to the depth of the rubber solid elements (1.67 mm). The latter does not present numerical difficulties since the Constrained_Lagrange_in_solid concept deals only with the node variables but nevertheless, the replacement of the necessary number of cord layers with fewer ones can change the overall inertia properties of the projectile and consequently its flexibility under dynamic loading. A reduction of the size of the solid elements certainly leads to more accurate results. The other disadvantage stems from the limited element dimensions that can be accommodated since the mesh arrangement is not arbitrary but must be related to the cord angles. The above limitations of the discrete model can present difficulties when used for large tyre projectiles. Having in mind the described features of the material model, a limited increase of the element size can be recommended. Once the size of the solid elements is decided a scaling coefficient for a coarse mesh, b, can be obtained for the cable material characteristics. For larger solid elements than the ones used in the model for the small projectile (while preserving the reinforcing angle of approximately 181): b ¼ ðd y =f z Þcoarse mesh =ðd y =f z Þfine model mesh 40

(6)

in order to maintain the reinforcement fraction characteristic of the coarse mesh the same as for the refined material model. In Eq. (6), d y is the distance between the reinforcing cords in the y-direction and f z is the number of cord layers per unit thickness (z-direction) associated with the ‘coarse mesh’ and the refined mesh of the material model proposed in this paper. For the calibrated material model

Fig. 19. Finite element mesh for full-scale tyre fragment simulation. Cords are given reduced diameter for clarity.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

665

ðd y =f z Þfine model mesh ¼ 9:365 mm2 . The cord diameter is then scaled as the square root of b. This approach was used to model a full scale, industry standard, tyre fragment test [24]. A tyre fragment, of dimensions 425  100  28 mm3 and mass 1.335 kg (wasted in the central section) is folded and placed in a large gas gun. When the fragment exits the gun barrel, it un-wraps and impacts the aircraft structure at a velocity of

Fig. 20. Comparison between experiment [25] and simulation [26] for full-scale wing access panel test: impact velocity 110 m s1 and impact angle 301.

ARTICLE IN PRESS 666

D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

114 m s1 and in the orientation required. The aircraft structure of interest here was a wing access panel [24,25]. Fig. 19 shows the tyre fragment mesh, where each element is 7.46  7.2  3.5 mm3, which should be compared to the previous case of 2.5  2.5  1.67 mm3 (Fig. 11). For the selected coarse mesh, the reinforcement angle is 17.71 and ðd y =f z Þcoarse mesh ¼ 28:572 mm2 taking into account that seven layers of cords are placed inside the projectile. Thus, b ¼ 3:047 for the current discretisation. In order to maintain the strength characteristic of the reinforcing layer unchanged, the cross-sectional area of the cords is scaled by coefficient b , which results in a cord diameter of 3.4 mm while the elastic modulus for the cable material remains the same. Fig. 20 compares the simulation [26] with experimental test [25] for an impact angle of 301. It can be seen that there is good comparison in respect of gross deformation fields. Ref. [26] also discusses strains at various positions on the wing access cover, and shows good agreement between experiment and simulation. It should be noted that the full-scale model was demanding computationally, taking over a day on a large Sun server. 9. Concluding remarks A numerical model for tyre fragment impact has been developed using standard material models in DYNA and using an approximation for the fibre reinforcement. The tyre model is approximate given the necessity for computationally robust and efficient modelling for industrially relevant problems. The model has been validated for small scale tyre fragment tests and for full-scale tyre fragment tests. Therefore, it can be concluded that the proposed simulation approach is accurate for a wide variety tyre impact scenarios, although demanding in computing time. It should be noted that there are limitations to the Mooney–Rivlin rubber model used here, and that there are more accurate rubber models in DYNA, which, however require a larger number of constants. Also, there are more sophisticated elastomeric models in the literature, see for example [27]. In addition, it is not straightforward to obtain the stresses within the tyre fragment directly using the output variables in LS-DYNA and the present simulation says nothing about the possible debonding between the rubber matrix and the reinforcing fibres. In fact, the diameter of the cords is comparable to the depth of the rubber elements and although approximates correctly the anisotropic material properties, it does not even closely approximate the physical connection between the cords and the rubber in the actual material. If this was required, then individual fibres would have to be modelled or some equivalent damage theory developed (e.g. see [28]). However, all this sophistication would increase computational demand and may lead to numerical instabilities, e.g., strain softening. It is proposed that, in the short term, a more productive avenue would be to simplify the current model by increasing mesh size or using simpler finite elements, to reduce computational time for real problems. As far as the contact algorithm is concerned, the standard contact model in DYNA gives good quality answers and tyre fragment behaviour seems to be insensitive to friction conditions. An adaptive mesh for the plate would seem to be unnecessary. Acknowledgement The support through EU Framework 5, DG XII – Competitive and sustainable growth—Key Action 4— New perspectives in Aeronautics: Crashworthiness under high velocity impact (CRAHVI) Contract no. G4RD-CT-2000-00395, is greatly acknowledged. Thanks are also due to Martin Kracht and Christian Bergler, of Cadfem Gmbh, for help in the full scale simulations and Jean Philippe Gallard, of CEAT, for full scale testing data. References [1] Johnson W, Sengupta AK, Ghosh SK. High velocity oblique and ricochet mainly of long rod projectiles: an overview. Int J Mech Sci 1982;24:425–36. [2] Johnson W, Sengupta AK, Ghosh SK. Plasticine modelled high velocity oblique impact and ricochet of long rods. Int J Mech Sci 1982;24:437–55.

ARTICLE IN PRESS D. Karagiozova, R.A.W. Mines / International Journal of Impact Engineering 34 (2007) 647–667

667

[3] Zukas JA, Segletes SB. Numerical modelling of hypervelocity impact phenomenon with desktop computer. Adv Eng Software 1993;14:77–84. [4] Zukas JA, Gaskill B. Ricochet of deforming projectiles from deforming plates. Int J Impact Eng 1996;18:601–10. [5] Knight NF, Jaunki N, Lawson RE, Ambur DR. Penetration simulation for uncontained engine debris impact on fuselage-like panels using LS-DYNA. Finite Elem Anal Des 2000;36:99–133. [6] Mines RAW, McKown S, Birch RS. Impact of aircraft rubber tyre fragments on aluminium alloy plates: I experimental. Int J Impact Eng, in press, doi:10.1016/j.ijimpeng.2006.02.005. [7] Watanabe Y, Kaldjian MJ. Modelling and analysis of bias-ply motorcycle tires. Comp Struct 1983;17:653–8. [8] Noor AK, Tanner JA. Advances and trends in the development of computational models for tires. Comp Struct 1985;20:533. [9] Tabaddor F, Stafford JR. Some aspects of the rubber composite finite element analysis. Comp Struct 1985;21:327–39. [10] Ngo D, Scordelis AC. Finite element analysis of reinforced concrete beam. J Am Concrete Inst 1967;64:152–63. [11] Reese S, Raible T, Wriggers P. Finite element modelling of orthotropic material behaviour in pneumatic members. Int J Solids Struct 2001;38:9525–44. [12] Phillips DV, Zienkiewicz OC. Finite element non-linear analysis of concrete structures. Proc Inst Civ Eng 1976;61(Part 2):59–88. [13] Sussman T, Bathe KJ. A finite element formulation for non-linear incompressible elastic and inelastic analysis. Comp Struct 1987;26:357–409. [14] Helnwein P, Liu CH, Meschke G, Mang HA. A new 3D finite element model for cord-reinforced rubber composites application to analysis of automobile tires. Finite Elem Anal Des 1993;14:1–16. [15] Liu CH, Hofstetter G, Mang HA. 3D element analysis of rubber-like materials at finite strains. Eng Comput 1994;11:111–28. [16] Ranjbaran A. Mathematical formulation of embedded reinforcement in 3D brick elements. Commun NumerMethods Eng 1996;12:897–903. [17] Ruter M, Stein E. Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity. Comput Methods Appl Mech Eng 2000;190:519–41. [18] Gerhard A, Gasser TC. A viscoelastic model for fibre-reinforced materials at finite strains: continuum basis. Comput Aspects Appl 2001;190:4379–403. [19] LS-DYNA users manual. Livermore Software Technology Corporation, March 1998. [20] Treloar LRG. The physics of rubber elasticity. Oxford: Clarendon Press; 1975. [21] Green AE, Adkins JE. Large strain elastic deformation and non linear continuum mechanics. Oxford: Clarendon Press; 1960. [22] Daley JR, Mays S. The complexity of material modelling in the design optimisation of elastomeric seals. In: Boast D, Coveney VA, editors. Finite element analysis of elastomers. London: Professional Engineering Publishing; 1999. p. 119–28. [23] Wriggers P. Computational contact mechanics. Chichester: Wiley; 2002. [24] Certification testing of 2024. Aluminium underwing access panel – tyre burst impact. Structures and mechanical systems laboratory report 5056, British Aerospace Ltd. 2000. [25] Gallard JP. Report on soft debris impact tests (CRAHVI deliverable D542). CEAT Int. Ref. S-01/5976000/P6; 2003. [26] Birch RS, Bergler C, Kracht M, Karagiozova D, Mines RAW. Post test simulation of airliner wing access panel subject to tyre debris impact. In: fifth European LS-DYNA users conference, 2005. [27] Criscone JC, Douglas AS, Hunter WC. Physically based strain invariant set for materials exhibiting transversely isotropic behaviour. J Mech Phys Sol 2001;49:871–97. [28] Allix O, Hild F. Continuum damage mechanics of materials and structures. Oxford: Elsevier; 2002.