MATHEMATICAL MODELS OF WAVE PROPAGATION AND FRACTURE IN UNIDIRECTIONAL COMPOSITES*
G. Osharovicha, M. Ayzenberg-Stepanenkob and V. D. Kubenkoc
aBar-Ilan University,Ramat-Gan, Israel
bBen-GurionUniversity of the Negev, Beer-Sheva, Israel
cInstitute of Mechanics, National Academy of Sciencesof Ukraine, Kiev, Ukraine
Abstract. Mathematical models and computer algorithms are developed to calculate dynamic stress concentration and fracture in a unidirectional composite.A plane problem is considered. The composite consists of a regular system alternating extensible (fiber)and pliable (adhesive) layers. Acalculation technique is designed to precise computingdynamicstress discontinuities that appear under fracture. Results of fracture simulations are presented in the composite pre-stretched along the fibers.
1 INTRODUCTION
Unidirectional layeredcomposite plates having a wide range of applications may berequired to endure intense dynamic loadings ([1-5]). A set of dynamic problem have been studied on the basis of the homogenization approach averaging local properties of microstructure (see e.g. [5-8]). Within such approach, features inherent in impact processes are, as a rule, explored by methods and models of solid dynamics. Those methods are well developed for analysis of a longwave spectrum, for which microstructure peculiarities, however, play a secondary role. Wave front gaps and high-gradient components responsible for dynamic stress concentrations and brittle fracturepropagation in a composite have not yet been adequately examined and remain a subject of contemporary research (see e.g. [9-10]). Significant rise in influence of microstructure on the wave picture essentially restricts capabilities of analytical modeling.There exist a set of numerical approaches intended for comprehensive simulations of composite dynamics in practice [11-14]. The recent work [14] reviews new software upgrades that are extending the digital advantage in the composites marketplace.
At the same time, explicit-time algorithms used in computer codes come across specific obstacles, which do not allow to calculate accurately wave fronts and high-gradient components. We mean the mesh dispersion responsible forthe emergence of short-wave "parasite" oscillations at calculation of transient processes and, notably, of fracture dynamics in reinforced composites: breaking of fibers and cracking of adhesive result in appearing of multiple wave fronts and high-frequency oscillations.
In this work, computational algorithms are designed on the basis of the mesh dispersion minimization (MDM) technique, which allowsmesh dispersion drawbacks to be eliminated or significantly decreased. The idea behind MDM is to properly adjust the so-called domains of influence determined by continuous and discrete models. The MDM approach proposed in [15] and upgraded in [16] has recently been used in diverse dynamic problems: high-speed penetration of layered shields [17], impact indentation of a rigid body into an elastic medium [18], and resonant excitation of lattice structures [19].
______*This research was supported by The Israel Science Foundation, Grant No. 504/08, and the Center of Advanced Studies of Mathematics at Ben-GurionUniversity of the Negev.
2 MECHANICAL AND MATHEMATICAL MODELS
2.1 Physical Assumptions. Description of Fracture Events
We consider elastic waves and dynamic fracture in a platemodeled by a periodic structure shown in Fig. 1.
The elastic plate of a unit width consists of an unbounded periodic system of high strength fibers, which alternate with pliable adhesive layers. Parameters of fibers: width, h, Young modulus, E, and density,; the adhesive has width, H, shear modulus, G, and density,. Axis x is directed across the fibers, and axis y – along them. In infinity (), fibers are stretched along the axis y by constant tensile stress .
It is assumed in the mechanical model that fibers function in tension-compression, while the adhesive is under shear stress (only longitudinal displacements exist in the plate). The assumption that normal stresses exist only in fibers while tangential stresses only in the adhesive is often used in studying static equilibrium of unidirectional composites. Such structures have wide range of practical application (for example, in aircraft and ship engineering). Although the stress state of structure components is, in fact, more complex, such approach correctly expresses the concept of the efficient performance of a reinforced material: high strength fibers are oriented along the tensile stress lines, while the adhesive facilitates a more uniform distribution of these loads between fibers, preventing stress concentrations.
The following problem is considered.At time all fibers are stretched along the axis y by a constant tensile stress applied in infinity (). Let one of fibers (say fiber) be suddenly fractured ina cross section (say in) at zero moment of time, t = 0see Fig. 1(b).As a result, the initial split increases with time Fig. 1(c). The problem is to describe the dynamic concentration of stresses, subsequent fracture process occurring with time and to reveal the “trauma” area after the (possible) fracture arrest.
We assume hH and introduce local coordinate X varied in interval , then , . Displacements and stresses in fibers are denoted by and , in adhesive and . The initial stressed state of the fibers is (). It was also assumed that initial tensile stresses in fibers is less than given critical value: . If , the initial fracture can excite a consequence breaking of fibers. Besides, we introduce the critical value of shear stresses in adhesive, : if , the shear crack appears in the corresponding point of adhesive and it can propagate away from the area of the initial fracture. At , the physical pattern of the process can be developed due to the following scenario:
at t = 0, two free boundaries appear in the fractured cross-section: and . At t > 0, the unloading wave begins propagating in the fiber, and its free ends move along axis y in opposite directions;
along with unloading, the broken fiber pulls adhesive adjacent to this area. Under the condition that adhesive does not reveal resistance to normal stresses, the line y = 0 in adhesive begins to open free: a transverse crack appears in adhesive between fibers m = 1 and m = 1, and the splitting area increases with time resulting in the intensification of shear wavesin adhesive, see Fig. 1(b);
the shear waves reach fibers nearest to broken fiber m = 0, overload them (i.e. in addition to ) and involve them in motion. In their turn, fibers excites shear waves in adhesive that propagate upward, to the next fibers, and backward, to the broken one. Shear waves reflected from the first intact fibers reach the broken one and decelerate it, while next intact fibers begin to overload. Together with this, free ends in fibers and shear cracks in the nearest adhesive can be appeared if normal and shear stresses are to be exceeded critical values * and *, respectively. These boundaries surround a domain associated to the trauma despite the fact that in some cases such trauma is not continued but can be formed from a set of alternative fractured and intact domains.
the dynamic process is developed with time farther on according to the described scheme, more and more volume of fibers neighboring to broken one are being fractured and involved into the motion, while intact those get additional load. At some particular moment of time, energy of the initial fracture is completely spent on tension-compression waves in fibers and shear waves in adhesive propagated away of the initial impact area. After that, strains and displacements of the composite reach a new static state, and dynamic process is finished. As a result, a trauma remains in this new static state.
The used model of the fiber dynamics describes the one-dimensional wave process in a thin rod embedded into adhesive, while two models are considered for adhesive:
(i)Model 1. This model is simplified one. It is designed under the assumption that adhesive can be represented by inertionless bonds of constant stiffness K=G/H. Note, the inertia of bond can be taken into account due to continuous distribution of the bond mass along the fiber in each the fiber-adhesive layer.
(ii)Model 2. This model is more precise. The adhesive is described by inertial bonds perceived shear stresses, while tension-compression stresses in bonds are neglected. Such a theoretical treatment of the nature of the components performance is partially justified by the fact that the shear modulus of adhesive is hundreds times less than that of fiber (see e.g. [1-4]), while their stretches have roughly the same level due to the cohesion of the fibers and the adhesive.
In Model 1, motion of the composite is completely described by displacements of fibers, while displacements of adhesive is calculated from the linear dependence on X: .
The following cases can be realized depending on strength conditions:
the composite remains intact: maximal stresses reached in fibers and adhesive do not exceed critical limits: . These conditions are used to calculate the parameters of stress concentration in intact fibers and adhesive associated to them;
fibers are fractured, while adhesive remains intact: ;
adhesive is fractured, while fibers (with ) remain intact: ;
both fibers and adhesive are fractured: .
2.2 Mathematical Formulation of the Problem
In the mathematical sense, the problem is related to a non-linear hyperbolic one possessing non-classical boundary conditions. Due to the symmetry, a quarter of plane x,y is considered in the calculation algorithm (let it be ). Displacements and strains of the composite at the intact static state (t < 0) are
,(1)
where is the strain in mth fiber. Fracture of 0th fiber at t = 0 changes (1) by adding condition :
. (2)
Let us reformulate the problem for the additional dynamic state. For this we subtract the static strains (1) from (2). Then boundary conditions for strains in fibers are the following:
(3)
The motion of fibers is described by 1D wave equation
. (4)
where correspond to reactive shear forces at the fiber-adhesive interface on the right and the left, respectively.Their expressions in Model 1 are:
, (5)
while in Model 2, where adhesive is considered as a array of inertial bond modeling adhesive, reactive forces are the following:
,(6)
and displacements in adhesive described by wave equations
(7)
with the following boundary conditions:
(8)
We add fracture conditions to system (1) (8) as follows. If in current cross section of fiber at the moment of time tension stress reaches critical value , this cross section breaks, and new boundaries appear resulting in the following conditionsfor equations (4):
(9)
In the case of the fiber-adhesive splitting, we, inaddition to (7), have the additional expressions for reactive forces in equations (4) and boundary conditions for :
,(10)
where + = 0, = H, and indexes “” at denote right and left interfaces, respectively.
3 MDM CALCULATION ALGORITHMS
To calculate the system (1) – (10) we use the explicit finite difference algorithm. Let temporal and spatial steps of the difference mesh be t, x and y. Let also values be measurement units (then is the velocity unit).
Calculation algorithms for a discrete analog of system (1) – (10) are constructed on the basis of the mesh dispersion minimization (MDM) technique elaborated in [15,16]. A simple example below allows the main principle of the MDM to be elucidated.
We mean the wave problem for semi-infinite straight elastic fiber subjected to the step tensile stress, which affects the free end. The equation of motion with the boundary condition, (I), the dispersion relation for free waves, (II), and the D’Alambert solution for the stress propagating along the fiber, (III), are:
, (11)
while H(z) is the Heaviside step function, c is the phase velocity. Initial conditions are zero. Relation (II) proves the dispersionless wave propagation: the phase velocity is independent on the wave length.
The discrete analog of equations (I) and (II) located in (11) is
.(12)
where q is the wave number, and indexes k and i are current numbers of temporal and spatial steps within the mesh . Note, we did not succeed in obtaining a closed analytical solution of Eq. (I) in (12) similar to (III) in (11).
In thediscrete model (12),c depends on q: waves propagate dispersionally. However, if = 1, the influence domains in the fiber,x =c0t, and its discrete analog (12) III, ix =c0kt, coincide. In other words, if = 1the mesh dispersion is eliminated. Indeed, in this case, the solution of Eq. (I) in (12) obtained by the mathematical induction is as follows:
(III) , (13)
whichcoincides with the D’Alambert solution in mesh nodes. Note, equality = 1 is the limiting value of the Courant stability criterion:.
In this simplest example, the MDM requirement, = 1, eliminates parasite effects of mesh discretization. It is a generalized concept of the well-known Courant condition that correlates dispersion properties of discrete and continuous models.
The second example is a single semi-infinite fiber upon an elastic foundation: only 0th fiber remains in system (4), Model 1, while : The fiber is loaded by rectangular pulse at its edge y= 0. Let , then . Stresses in the fiber shown in Fig. 2 (a) are calculated by algorithm (12) withxandtReactive forces are approximated by conventionaldifference form .In Fig. 2, (b) the same stressesare calculated by MDM algorithm, in which tx= 1, and the three-point approximation of reactive forces is made:
.(14)
Representation (14) allows domains of influence in continual and difference problemsto be coincided. Comparing results one can see spreading fronts and parasite oscillations saturating the wave picture in the former case. This fact prevents us from calculating the fracture problem on the basis of the stress criterionconsidered above. Note, stresses shown in Fig. 2 (b) coincide with those obtained from the analyticalsolution [20]: the MDM algorithm can be successively used for the correct description of fracture.
4 RESULTS OF COMPUTER SIMULATIONS
In Model 1, the MDM calculation algorithm for an array of fibers is proved if to set cftx in Eq. (4), while shear deformations in adhesive, Eq. (5), are approximated as
(15)
As to Model 2, the following approximation results in the MDM algorithm (t = 1):
Some results of computer simulations conducted on the basis of MDM algorithms (15) and (16)are presented below. Parameters of the fiber are measurement units: h= 1, E= 1, (cf = 1); adhesive parameters are H = 5, G = 0.025, a = 0.4. We set also , and ty = 1, xca within the MDM algorithm.
First, let us discuss Model 1. In Fig. 3, stresses vs. time are depicted for the case in which the composite remains intact at t > 0: andare less than corresponding maximal values reached after the initial fracture at t = 0. They are observed in cross-section y = 0: . Note, if t > 25H/cf, stresses are closed enough to their static limits: . Below, in calculation of fracture problems, intervals between maximal and static values serve for setting and.
Fracture development with time is shown in Fig. 3. One can observe processes of propagation and arrest of normal cracks – (a) and (b), and shear cracks – (c).
Results shown in Fig. 5 are calculated in the casewhen composite remains intact at t > 0. The gap of shear stresses, initially equal to 0.1, doubled after the first reflection from intact fiber m = 1 and remains equal to 0.2 at the whole process. After next reflections, duration of the gap decreases, and influence of gaps to normal stresses decreases, as well.As comparison shows, the maximal and static stresses in fibers are practically the same as those obtained for Model 1. But the significant difference in maximal shear stresses (~ 3 times) will result in different fracture processes. Thus,Model 1can be used at the initial, pre-fracture stage, where a stress concentration pattern is calculated, while fracture propagation processes saturated by front gaps are to be computed with Model 2.
It could be used within computer simulators allowing optimization problems to be explored. Let, for example, limit be constant, while limitbe variedwithin interval. Our aim is to find such a value of , at which the volume of breaking fibers,(or, that is the same, the normal crack length equal to 2H), is less than given value M. In the table below we present some related results obtained in the case and M = 10; in two lower rows we set the amount (its half due to symmetry) of fractured fibers, , and lengths of shear cracks in adhesive, . One can see that the required result is . Note, in presented example, shear cracks are observed only in layer .
/ 0.07 / 0.075 / 0.08 / 0.085 / 0.1 / 0.15 / 0.25/ 0 / 0 / 5 / 7 / 9 / 13 / 13
y*/h / / 208 / 65 / 23 / 11 / 3 / 0
Fractureprocess in the case possesses the time pattern shown in Fig. 6.
REFERENCES
1. Gibson R F, Wilson D G: ‘Dynamic mechanical properties of fiber-reinforced composite materials’. Shock and Vibration Digest 1979 11(10) 3-11.
2. Sela N, Ishai O: ‘Interlaminar fracture toughness and toughening of laminated composite materials: a review’. Composites 1989 20(5) 423-35.
3. Ashbee K: Fundamental principles of fiber reinforced composites. Technomic. Lancaster, PA, 1993.
4. Deng S, Ye L, Mai Y: ‘Influence of fibre cross-sectional aspect ratio on mechanical properties of glass fibre/epoxy composites. 2. Interlaminar fracture and impact behaviour’. Composites Science & Technology 1999 59 1725-34.
5. Achenbah J D: Vibrations and waves in directional composites. In: Mechanics of
Composite Materials 2. N-Y. Academic Press, 1975.
6. Achenbach J D, Sun C T, Herrman G:'Continuum theory for a laminated medium'. J. Appl. Mech. 1968 35 467–75.
7.Santosa F, Symes W: ‘A dispersive effective medium for wave propagation in periodic composites’.SIAM J. Appl. Math. 1991 51 (4) 984-1005.
8. Weekes S: ‘A stable scheme for the numerical computation of long wave propagation
in temporal laminates'. J. Comput. Phys. 2002 176 345–62.
- Dai J J, Liang X D, Yao X F, Yeh H Y: ‘Study of cracked unidirectional glass fiber-reinforced composites by digital speckle correlation method’. J. Reinforced Plastics and Comp. 2005 24 (16) 1737- 46.
- Nayak A K Shenoi R A, Moy S J: ’Transient response of initially stressed composite plates’. Finite Elements in Analysis and Design 2006 42 (10) 821-36.
11. Langlie S, Cheng W, Dilber I: ‘Computer simulation of high-speed impact response of composites'. Proc. 1990 ACM/IEEE Conf. on Supercomputing,N-Y, 1990.
12. Reddy J N, Robbins D H Jr: ‘Theories and computational models for composite laminates'. Appl. Mech. Rev. Trans. ASME 1994 47 147-69.
13. Silling S A, Taylor,P A: ‘The simulation of composite material response under dynamic compressive loading’. Modelling Simul. Mater. Sci. Eng. 1994 2 689-99.
- McConnell V P: ‘Software delivers faster time-to-part and reduced testing’. Reinforced Plastics 2007 51(8) 24-8.
15. Abdukadyrov S, Alexandrova N, Stepanenko M: ‘On MDM approach to calculation of solid and structure dynamics’. J. Mining Sci. 1984 6 34-40.
- Ayzenberg-Stepanenko M: ‘Wave propagation and fracture in elastic lattices and directional composites’. Personal Armour Systems - Conf. Proc., Colchester, British Crow Copyright/MOD, 1998.
17. Ayzenberg-Stepanenko M, Osharovich G: ‘Penetration models for optimization of composite armor shields under high-speed impact’. Sixth Russian-Israeli Bi-National Workshop. Jerusalem, IsraeliAcademy of Sciences and Humanities, 2007.
18. Kubenko V D, Ayzenberg-Stepanenko M V: ‘Impact indentation of a rigid body into elastic layer. Analytical and numerical approaches’. Math.Methods and Phys. Mech. Fields 2008 51(2) 61-74.
19. Ayzenberg-Stepanenko M, Slepyan L: ‘Resonant-frequency primitive waveforms and star waves in lattices’. Journal of Sound and Vibration 2008 313 812-21.
- Slepyan L: Transient elastic waves. St-Petersburg, Sudostroenie, 1972 (in Russian).
1-1