International Journal of Fracture 128: 325–333, 2004. © 2004 Kluwer Academic Publishers. Printed in the Netherlands.
Numerical simulation of brittle fracture of thin-walled structures1 A.V. SHUTOV Lavrentyev Institute of Hydrodynamics SB RAS, 630090 Novosibirsk, Russia
Abstract. The model of brittle material, allowing for the determination of initiation and growth of cracks in thinwalled structures up to their destruction, is offered. Solutions of some problems are submitted: destruction of a plate at a homogeneous stress state; destruction of a beam at a pure bending; dynamic deformation and destruction of a plate with the concentrated mass having various initial speeds. Destruction of part of the material of the plate is shown to result in dynamic loss of stability of the solution. Key words: Brittle fracture, dynamic buckling, finite elements, shell.
1. Introduction This work constitutes further inquiry into the nonlinear material model proposed earlier (Suidan and Schnobrich, 1973; Bathe and Ramaswamy, 1979). As brittle fracture goes, we understand fracture without taking into account viscosity and accumulation of microdamages. It is supposed that the criterion of fracture is the excess of the maximal main stress over a limit tensile strength of a material. The material at achievement of strength instantly passes from the elastic into the destroyed state. This model of brittle material is entered into the library of material models for a shell element (Korobeinikov and Bondarenko, 1996) of PIONER code (Korobeinikov et al., 1989). With the help of this code, solutions both of test problems and new problems on the destruction of bronze-plate protective armour of ancient nomads are obtained. 2. The numerical solution of nonlinear problems on the deformation of thin-walled structures We shall write down the statement of a problem and a derivation of the nonlinear discrete equations of shell theory used in this work. More detailed derivation of the discrete equations of shell theory is given, for example, in Bathe (1982), Zienkiewicz and Taylor (1991), Korobeinikov and Bondarenko (1996). 2.1. T HE EQUATIONS OF DEFORMATION OF SOLIDS Consider that deformations, rotations and displacements are small. We shall obtain the basic equations of the problem. 1 The work is supported by RFBR (codes of projects 02-01-00195, 04-06-80248) and by the grant of the President
of the Russian Federation (SS319.2003.1).
326 A.V. Shutov 2.1.1. The dynamic equilibrium equation with boundary and initial conditions ¨ ∇ · σ = ρ u,
u = u∗ at Su ,
N · σ = T∗ at ST ,
u|t =0 = u0 ,
˙ t =0 = v0 . (1) u|
Hereinafter σ Cauchy stress tensor; ∇ · σ stress tensor divergence; u displacements; V the body; S the boundary of V ; Su , St parts of a surface S = Su ∪ ST at which displacements u and the stress vector T ≡ N · σ are given respectively; N the unit external normal to ST ; u0 , v0 the given vector fields; ρ mass density of a material; the point above a value designates a partial derivative of this value in time t (deformation parameter). 2.1.2. Strain-displacement relations (ε is Cauchy strain tensor) 1 ε = (∇u + ∇uT ). 2 2.1.3. Constitutive law σ˙ = C : ε˙ ,
(2)
where C, the fourth rank tensor whose components generally depend on a history of deformation. The equations of quasistatic deformations are derived by neglecting the inertial members in the right part of Equation (1). We shall consider the weak form of the movement equations Equation (1) (σ : δε + ρ u¨ · δu)dV = T∗ · δudS ∀ δu, (3) V
ST
where δu is admissible virtual displacement field. The step-by-step procedure of integration in time of Equation (3) is used. Assume that at the moment t all required values are determined. The appropriate values at the moment t + t must be determined. Equation (3) at the moment t + t yields t +t t +t ∗ t (σ : δε + ρ T · δudS − σ : δεdV ∀ δu. (4) u¨ · δu)dV = V
ST
V
Hereinafter the symbol means increment of a value from the moment t up to the moment t + t. Linearizing relations in Equation (2) concerning the moment of time t, we obtain σ = t C : ε.
(5)
Substituting Equation (5) into Equation (4), we obtain the linearized motion equation t +t t t +t ∗ t (δε : C : ε + ρ T · δudS − σ : δεdV ∀ δu. u¨ · δu)dV = V
ST
(6)
V
After FEM discretization the ODE system (6) is solved step-by-step by Newmark’s method. As the system (6) is linearized the solution can be specified by iterative procedure. 2.2. F INITE ELEMENTS FOR CALCULATION OF THIN BODIES If one linear size of an element essentially surpasses another (others), the tangent stiffness matrix (or the effective stiffness matrix) becomes ill conditioned. This can be avoided by entering the static Kirchhoff-Love hypothesis: the stress normal to the midsurface of the shell is negligible. To introduce this hypothesis into curvilinear three-dimensional finite elements,
Numerical simulation of brittle fracture of thin-walled structures 327 the local coordinate system in an element is guided in such a manner that the coordinate t is directed along ‘normal’ to the midsurface. The word ‘normal’ is quoted, as only the affinity of this vector to the real normal is required. A layer is defined by a surface t = tint , where tint –the value of local coordinate t at an integration point. At each integration point two tangent vectors to a layer of unit length eˆ 1 , eˆ 2 , orthogonal to each other, are defined. The third vector of unit length eˆ 3 is entered as a normal vector to a layer, normal to the first two. We shall designate the local Cartesian coordinate system with orthonormal basic vectors eˆ 1 , eˆ 2 , eˆ 3 as the layer coordinate system. Introduction of the static Kirchhoff-Love hypothesis corresponds to the assumption σˆ 33 = 0. We enter vectors σˆ and εˆ : T T σˆ = σˆ 11 , σˆ 22 , σˆ 12 , σˆ 13 , σˆ 23 , εˆ = εˆ 11 , εˆ 22 , 2ˆε12 , 2ˆε13 , 2ˆε23 . Now we can rewrite Equation(5) by introduction of restriction σˆ 33 = 0 in the form ˆ ε, σˆ = t Cˆ ˆ is a matrix. Along the t direction at this element, linear interpolation is used. where t C The element considered has a zero strain energy spurious eigenmode. For suppression of this mode except for a static hypothesis we shall consider the updated kinematic Timoshenko hypothesis (a fibre directed along a normal to the midsurface of the shell is not bent and stretched during deformation). Updating this hypothesis consists of, instead of a fibre directed along a normal to the midsurface, our use of a fibre connecting two nodes, on the bottom (t = −1) and top (t = +1) surfaces of an element. It allows for simplified finite-element modeling. 3. Model of a brittle material of a thin body The behavior of a brittle material is determined by three constants: E, Young modulus; v, Poisson ratio; σt , tensile strength of a material. We introduce three parameters determined at each point of a body for each moment of time t: nic , parameter of i-th (i = 1, 2) crack (it is supposed that at each point of a body one crack or two cracks can exist simultaneously). nic = 0 if material is undamaged; nic = 1 if the crack is opened; nic = 2, if the crack is closed. Let ϕ, angle between the first stress principal direction and eˆ at the moment of crack occurrence. At the initial moment (t = 0) the material is supposed intact and at all integration points 0 i nc = 0 (i = 1, 2) is given, and the parameter ϕ is not defined. For the intact material, Hooke’s constitutive law with the matrix, which takes into account the static Kirchhoff-Love hypothesis, is used (k = 5/6). 1 v 0 0 0 1 0 0 0 1−v E 0 0 E . ˆ (7) C = 2 2 1−v 1−v sym. 0 k 2 1−v k 2
328 A.V. Shutov Let the values of parameters t nic , ϕ at time t and strains t +t εˆ at the moment t + t be ˆ and a vector t +t σ . known. We shall write out the algorithm of definition of a matrix t +t C t i Item 1. The material is not damaged ( nc = 0, i = 1, 2). Trial stresses are defined ˆ E t +t εˆ . σ˜ˆ = C Influence of σˆ 13 , σˆ 23 on formation of a crack is neglected (this assumption is carried out precisely with the use of the kinematic Kirchhoff hypothesis). For definition of an opportunity of occurrence of a crack (cracks) we find principal values and principal directions of a stress tensor σ˜ˆ at a layer plane and angle ϕ between the first principal direction and eˆ 1 . We shall designate main trial stresses at a layer plane through σ˜ 1 , σ˜ 2 , σ˜ 1 σ˜ 2 . Compare σ˜ 1 with σt . 1.1. If σ˜ 1 < σt , at the given point the material is supposed to be intact. Finally we determine t +t
ˆ =C ˆ E, C
t +t
σˆ = σ˜ˆ .
i
Suppose t +t nc = 0 (i = 1, 2), then the angle ϕ is not remembered. 1.2. If σ˜ 1 σt and σ˜ 2 < σt , it is supposed that at the given point, the crack, whose plane is normal both to the layer surface and to the first principal direction, is formed. The crack is modeled by orthotropic elastic material, weakened in the first principal direction. We update stresses σ1 = 0,
σ2 = σ˜ 2
σ12 = 0.
Stresses t +t σˆ at a layer coordinate system are obtained by transformation from stresses determined in the principal direction coordinate system t +t
σˆ 11 = σ2 sin2 ϕ,
t +t
σˆ 13 = k
t +t
E t +t εˆ 13 , 1+v
σˆ 22 = σ2 cos2 ϕ, t +t
σˆ 23 = k
t +t
σˆ 12 = −σ2 sin ϕ cos ϕ,
E t +t εˆ 23 . 1+v
ˇ is defined as For orthotropic material, in modeling a crack, the matrix C 0 0 0 k1 vk1 1 0 0 0 1−v E k 0 0 2 . ˇ = C 2 2 1−v 1−v 0 sym. k 2 1−v k 2
(8)
ˆ is obtained by rotation from a matrix C. ˇ Constants k1 , k2 in Equation (8) The matrix t +t C (small positive numbers) are required for regularization of a tangent stiffness matrix t +t K. 1 Suppose t +t nc = 1, and we remember the value of ϕ which further does not change. 1.3. If σ˜ 1 σt and σ˜ 2 σt , material at the given point is completely destroyed, so
t +t
σˆ = 0,
Numerical simulation of brittle fracture of thin-walled structures 329 0 0 0 k1 vk1 0 0 0 k1 1−v E k 0 0 2 . ˇ = C (9) 2 2 1−v 1−v k2 0 sym. k 2 1−v k2 k 2
ˆ is obtained by rotation from a matrix C. ˇ Suppose t +t nic = 1, (i = 1, 2) the Matrix t +t C value of the angle ϕ is not important. Item 2. The material has only one crack (t n1c = 0, t n2c = 0). We determine t +t ε in the coordinate system connected to orientation of a crack. We find trial σ˜ 2 from Hooke’s law σ˜ 2 =
E (vε11 + ε22 ). 1 − v2
Check an opportunity of occurrence of a crack, whose plane is orthogonal to the plane of the first crack, comparing value σ˜ 2 with strength σt . 2.1. If σ˜ 2 σt , at the given point there is a second crack normal to the first crack, and i Equation (9) is valid. Suppose t +t nc = 1, (i = 1, 2), then the value of the angle ϕ is not important. 2.2. If σ˜ 2 < σt , at a point there is one crack. We examine a crack: whether the crack is open or it is closed. ˆ are 2.2.1. If ε11 < 0, it is supposed that the crack is closed. Stresses t +t σˆ and matrix t +t C defined from Hooke’s law t +t
ˆ E t +t εˆ , σˆ = C
t +t
1
2
ˆ =C ˆ E. C
Suppose t +t nc = 2, t +t nc = 0, then the angle ϕ does not change. 2.2.2. If ε11 0, it is supposed that the crack is open. The finding of stresses t +t σˆ and ˆ is carried out just as in item 1.2; suppose t +t n1c = 1, t +t n2c = 0, then ϕ does matrix t +t C not change. Item 3. The material is destroyed with loss of carrying capacity (t nic = 0, i = 1, 2). It is supposed that at full destruction the material already under any conditions can not restore ˆ are defined just as in item 1.3; suppose carrying capacity. Stresses t +t σˆ and matrix t +t C t +t i nc = 1(i = 1, 2), then a value ϕ is not important. 4. Results of numerical calculations The algorithm given in Section 2 is realized in the model of brittle material entered into the library of models of materials for finite element of a shell of the computer code PIONER. 4.1. P LATE AT PLANE STRESS We shall model the destruction of a plate at stretching in two orthogonal directions in quasistatics (Figure 1).
330 A.V. Shutov
Figure 1. Plate at plane stress: system configuration and boundary conditions (a); prescribed displacements (b); resulting stresses (c).
Figure 2. Destruction of a beam at a pure bending: system configuration and boundary conditions (a); fracture zone expansion (b); distribution of σxx stress (c).
The plate is stretched in two orthogonal directions by the given displacements u and u/2 at its edges (Figure 1a). The following constants are used: E = 103 GPa, v = 0.25 GPa, σt = 200 MPa. Boundary conditions are put so that in the plate the plane stress is carried out. The given displacement u grows linearly from zero up to umax (Figure 1b). The solution of a problem is submitted in Figure 1c. 4.2. D ESTRUCTION OF A BEAM AT A PURE BENDING Consider the problem of destruction of a beam with rectangular cross section under pure bending in quasistatics. Such bending can be carried out both by the moment M, and by the angle of a rotation β of an edge around of the axis (Figure 2a). Under pure bending of an elastic beam the analytical solution of this problem is known. The σxx and bending moment M can be defined through angle β under the formulas EI Ez β, M= β. (10) L L Here E Young modulus, L half of length of a beam, I the moment of inertia of cross section of a beam. From Equation (10) we find critical values of angle βcr and the appropriate moment Mcr , at which on a surface z = h/2 (h, the height of a beam) destruction of a material begins (σxx |z=h/2 = σt ) σxx =
Numerical simulation of brittle fracture of thin-walled structures 331
Figure 3. Numerical solution on beam destruction: FEM mesh (a); bending moment M versus bending angle β (b,c,d).
βcr = 2Lσt /(Eh),
Mcr = 2I σt / h.
(11)
While β βcr , the solution (10) is valid. At β > βcr part of the material is destroyed (it is shaded on Figure 2b), and the rest (an elastic core) continues to deform elastically. Distribution of σxx stress is shown in Figure 2c. Moment M for the given angle β βcr is the following M=
8I σt3 L2 1 . h3 E 2 β 2
(12)
For the numerical solution we consider half of beam. At the left edge are symmetry conditions and at the right edge the corner β of a bending is set. Geometrical parameters of a beam are given on Figure 3a. The beam is modeled by one shell element with five degrees of freedom at nodes with cubic interpolation along an axis of a beam and linear interpolation along a cross direction. If the full order of numerical integration 4 × 2 × 2 on coordinates r, s, t is used, the analytical solution (10) is reproduced precisely. However, two integration points along the thickness of a beam will be insufficient for the solution of a problem in case of destruction of a material. We use the following constants: E = 300 GPa, v = 0, σt = 100 MPa. According to Equation (11), critical values βcr = 0, 008 and Mcr = 1667 N cm are obtained. The dependence of the moment M on angle β, obtained in the exact solution according to Equation (12), is shown in Figures 3b–d by a continuous line. Numerical solutions with the use of different numbers of integration points along coordinate t are submitted by points. It is apparent that 16 integration points along the thickness of a beam already reproduce the analytical solution well enough. 4.3. DYNAMIC DEFORMATION AND DESTRUCTION OF A BRONZE PLATE WITH THE CONCENTRATED MASS ADDED
Consider a problem on the destruction of a bronze plate, two edges of which are pinned, and other two are free (Figure 4). In the centre of the plate the initial speed v0 , directed on the normal to the plane of the plate, is given; everywhere else initial speed is equal to zero. Initial displacements are zero. At the centre of the plate (in that point where initial speed is given), the concentrated mass m = 0.08 kg, modeling mass of arrow, is entered. Geometrical parameters of the plate are
332 A.V. Shutov
Figure 4. System configuration, boundary conditions, and FEM mesh.
Figure 5. Numerical solution on plate destruction: central deflection w versus time t (a); deformed configuration and fracture zone (b).
Figure 6. Chaotic dispersion of plate fragments: (a) at time t = 1.42 × 10−4 s and (b) at time t = 2.54 × 10−4 s.
given in Figure 4. They are close to the sizes of plates out of which the armour of Central Asian nomads were made. The basic purpose of the solution of this problem is to reveal a character of destruction of the plate at various values of initial speed v0 . A material of the plate – bronze with: E = 103 GPa, v = 0.25, σt = 200 MPa, ρ = 8.8 × 10−5 N s2 /cm4 . The plate is modeled by nine-nodal finite elements. The solution is carried out with 3 × 3 × 4 numerical integration along coordinates r, s, t. We shall produce solutions of this problem for different values of initial speed v0 . Dependence of central deflection w from time and the deformed configuration of the plate, received in calculations with initial speed v0 = 3, 2 m/s, is given in Figure 5 (r is the scaling factor of displacements). On each element the number of integration points at which there were cracks (maximum 36) is given. At an increase of initial speed up to v0 = 10 m/s, the plate is completely destroyed. At destruction of a plate, loss of stability of the solution is revealed. According to the stated boundary conditions, the solution of the problem should have some symmetry.
Numerical simulation of brittle fracture of thin-walled structures 333 If the solution is stable in relation to small perturbations of coordinates of nodal points (in our case the perturbations are brought during automatic mesh generation), the appropriate components of translations should differ for all time of the solution a little. However, continuous dependence of the solution on perturbations of coordinates of nodal points is lost. For large enough times (Figure 6a corresponds to time t = 1.42 × 10−4 s, Figure 6b – to time t = 2.54 × 10−4 s), this asymmetry is obtained in the numerical solution. 5. Conclusion The nonlinear motion equations of thin bodies with FEM discretization are considered. The constitutive law of brittle material concerning material nonlinearity is offered. The model of brittle material allows for taking into account the occurrence of cracks, their closing or development up to the exhaustion of the structure’s bearing ability. This model is entered into the library of models of materials of finite elements of a shell of the PIONER computer code. Solutions of some test problems on the destruction of the plate, providing an opportunity to show the computer complex PIONER under the solution of problems on nonlinear fracture mechanics, are produced. The dynamic problem on the destruction of a bronze plate with the concentrated mass added to the given initial speed is solved. It is established that destruction of a plate is accompanied by the phenomenon of dynamic buckling. References Bathe, K.-J. (1982). Finite Element Procedures in Engineering Analysis. Englewood Cliffs, New Jersey, Prentice Hall. Bathe, K.-J. and Ramaswamy, S. (1979). On three-dimensional nonlinear analysis of concrete structures. Nuclear Engineering and Design 52(3), 385–409. Korobeinikov, S.N. and Bondarenko, M.I. (1996). A material and geometric nonlinear analysis of shells including large rotation increments. Numerical Methods in Engineering’ 96: Proc of the 2nd ECCOMAS Conf. (Edited by J.-A. Désidéri et al.). Chichester, John Wiley & Sons, 754–762. Korobeinikov, S.N., Agapov, V.P., Bondarenko, M.I. and Soldatkin, A.N. (1989). The general purpose nonlinear finite element structural analysis program PIONER, Proc. of the Intern. Conf. on Numerical Methods and Applications (Edited by B. Sendov et al.). Sofia: Publ. House of the Bulgarian Acad. of Sci., 228–233. Suidan, M. and Schnobrich, W. (1973). Finite Element Analysis of reinforced concrete. Journal of the Structural Division (Proceedings of the ASME) 99(10), 2109–2122. Zienkiewicz, O.C. and Taylor, R.L. (1991). The Finite Element Method. London et al.: McGraw Hill.