International Journal of Fracture, Vol. 13, No. 4, August 1977 Noordhoff International Publishing- Alphen aan den Rijn Printed in The Netherlands

443

Analysis of fast fracture and crack arrest by finite differences M. S H M U E L Y Technion -Israel Institute of Technology, Haifa, Israel

(Received October 5, 1976) ABSTRACT A finite difference scheme, recently proposed and mainly concerned with the DCB static solution under plane strain conditions, is here further extended to provide more information on fast crack propagation and arrest. The analysis is shown to predict crack lengths at arrest and crack velocities as well as dynamic stress intensity factors which are very close to those reported in the literature. This good agreement with experimental findings was obtained for steel and epoxy resin which, while in different classes of material, nevertheless as far as the analysis is concerned are distinguished only by their Poisson's ratios. 1. Introduction In [1], the finite difference method was shown to be applicable in solving elast o d y n a m i c p r o b l e m s of fracture mechanics. In particular, the solution of the propagating central crack had been found to yield results which are in good agreement with reality, both qualitatively and quantitatively. In [2], the method was further extended to include finite geometries and uneven load patterns, with the double cantilever b e a m specimen (DCB) serving as a test case. There, while being mainly concerned with the finite difference static solution, preliminary results for the dynamic case have been found to be in good a g r e e m e n t with experimental findings. H e r e , the dynamic solution corresponding to the DCB problem is further elaborated to obtain a clearer insight into the effects of the dynamic stress and energy distribution during fast crack propagation and at arrest. Available two-dimensional dynamic analyses in fracture mechanics [3, 4, 5] are too specialized to be applicable to the DCB problem. A b e a m on elastic foundation model recently p r o p o s e d [6] for handling the DCB dynamic problem was found to yield results which are in good agreement with experimental findings. Being a one-dimensional solution, h o w e v e r , with stress levels and displacements transversely averaged, the model can only a p p r o x i m a t e the actual situation, as stress w a v e reflections f r o m the boundaries and stress w a v e diffractions at the crack tip seem to play a significant role. The finite difference dynamic solution presented here, being spatially two-dimensional, readily responds to the last mentioned peculiar occurrences. As in [1] and [2], a fracture stress criterion is assumed, according to which the crack extends by one mesh size w h e n e v e r the cleavage stress at a given short distance f r o m the tip reaches a predetermined value. H e r e , as in [2], starting with exactly the same numerical initial conditions, solutions are obtained for different critical stresses which are taken to be smaller than or equal to the stress level ahead of the crack tip at initiation and are kept constant during each run. Here, in addition, we show that the critical stresses, being proportional to and representing the stress intensity factors, m a y be considered as the averages of those critical stresses which actually exist in each case during the propagation stage. The solution, then, predicts apparent dynamic stress intensity factors which are in good a g r e e m e n t with those a p p r o x i m a t e d or measured by other Int. Journ. of Fracture, 13 (1977) 443-454

444

M. Shmuely

means, when referring to the same crack velocities or, alternatively, given a particular initial geometry, to the same crack lengths at arrest. With the objective of comparing the numerical results with experimental findings, the analysis here is related to the two materials, epoxy resin and steel, which were separately employed in two experimental studies of the DCB conducted respectively at the Institut fiir FestkSrpermechanik (Freiburg) [7] and the Battelle (Columbus) Laboratories [8]. The same geometry was used in both studies and accordingly was introduced into the present analysis. In the numerical analysis employing dimensionless physical quantities, the two investigated materials (which are assumed to be linear elastic) are distinguished by the Poisson ratio only. It was interesting to find that although the numerical schemes used were identical except for the Poisson ratio, the final results, corresponding respectively to the two different materials, diverge so as to yield solutions which are in good agreement, each, with its experimental counterpart.

2. Elastic equations and numerical scheme

Only an outline of the equations and numerical procedures used is presented here. A detailed description of the finite difference approximations employed and the structure of the final numerical scheme is given in [2]. Assuming an elastic medium under plane strain conditions, the governing equations are:

(C~- C~)Uk.k, + C~U,.kk = U,.,, + aU,., i, k = l, 2

(1)

In Eqns. (1), C~ and C2 are respectively the dilatational and distortional wave velocities of the material considered, Ut stands for the displacement vector components in the ! direction, t is the time and the ordinary tensor notation for differentiation and summation is used. Depending on the value substituted for a, the same numerical scheme may be used for either the static case or the dynamic case. By letting a be equal to zero, Eqns. (1) reduce to the well-known elastic equations of motion, yielding the dynamic solution; whereas if a > 0, successive changes in displacements, occurring in real time, correspond to successive numerical two-point semi-iterative cycles, converging finally to the time independent version of Eqns. (1), i.e., the static solution. Hence, to simulate a complete fracture process, a single numerical scheme is employed first to yield the static solution (a = 0) and then to handle the fast crack propagation and arrest. Thanks to the symmetrical loading of the DCB (Fig. la) it suffices to consider only the upper part (y ~>0) of the specimen, provided that the boundary conditions related to that part, i.e.: tryy=crxy=0 tr~=o'~y=0

on on

x=O,L+e

y=0

for0~
(2a) (2b)

tryr=trxr=0

on

y=H

(2c)

and for0~
are supplemented by the symmetry conditions Uy=~rxy=0

on

y=0

for a o + e < x < - L + e

Int. Journ. of Fracture, 13 (1977) 443-454

(2d)

Analysis of fast fracture and crack arrest by finite differences

445

~Y

~(Xo;yo)

L

"t

---~X

tp

oo-

a(O~ u

(a) P

(b)

sz

(c)

s2

Figure 1. The DCB specimen configuration; (a) actual dimensions; (b)-(c) contour of loaded specimens as obtained from analysis.

Since, in practice, the load pins are essentially fixed during the propagating stage, a constant deflection 8 should be assigned to the point of load application (x0, Y0),yielding the additional condition

Uy(xo, Y0) = 8

(3)

Equations (1) together with conditions (2) and (3) define the problem completely, provided that, in the analysis, the constant deflection is assigned to a point which, while being kept at a horizontal distance of a0 from the crack tip, is located on the boundary of the specimen (Figs. lb and lc) rather than inside it. If, instead, the loaded point is moved to its actual position inside the specimen, a cut, which passes through that point, should be introduced into the otherwise simply connected region, and the appropriate conditions on the boundary of that cut should be added to the above mentioned system of equations. Employing the just described extended equation set would,only then be of any advantage if the cut is taken along the pin hole periphery, to obtain a solution closer to reality. For the conditions on the boundary of the relatively small hole to be properly approximated, a denser mesh than that which would have otherwise been used to obtain satisfactory results is required. It was found, however, that by reducing the mesh size, while in any other respect the solution is only slightly improved [2], the crack end shape is significantly changed, yielding steeper slopes closer to the analytical elliptic profile (see profiles $2 in Figs. lb and lc) for smaller mesh sizes. It was shown in [1] and verified here again that the dynamic solution is only then close to reality when the crack end shape is of the cusp type. While in [1], considering the homogeneously loaded central crack model, the crack end profile sought was imposed Int. Journ. of Fracture, 13 (1977) 443-454

446

M. Shmuely

by constraining the displacements at the crack tip to follow a prescribed pattern, here, thanks to the special loading configuration, the profile sought is self-obtained when the mesh is sufficiently large (see profiles Sl, Figs. lb and lc). It was consequently decided to tackle the problem with the first suggested equation system (i.e., the constant deflection assigned to a boundary point), while using the largest possible mesh size. The accuracy of the numerical analysis is checked via the static solution. It was shown in [2] and verified here again, that the stress intensity factor (K~) levels, as obtained numerically by means of the J integral [4], are close to those reported in the literature. A good matching was also found between the stiffnesses 8/P as predicted by the analysis and those measured in experiments. In [2], where the load was applied as illustrated in Fig. lc however, the best matching was obtained when the load pins in the experiments were small compared with those ordinarily used in DCB tests. Here, in addition, we found that a good match is also obtained with ordinary sized pins if the load in the scheme is applied as illustrated in Fig. lb. The numerical static solution was also compared with results obtained from a photo-elastic analysis. To this end, given the displacement field and the appropriate constitutive equations, the fields of normal, shear, principal and difference of principal stresses were calculated and recorded on disk. Given each of the last mentioned fields, a special program was used to trace through the field isostress lines. These lines are then written on the CALCOMP plotter. A detailed description of numerical procedures and experimental devices used, as well as results obtained, will be reported in the near future elsewhere. Two comparative examples of a more qualitative nature are shown in Figs. 2 and 3. In order to decide upon the largest possible mesh size, all lengths in the scheme are scaled in terms of the specimen height H (Fig. la). Results are first obtained for a

Figure 2. The analytic shear stress field compared with isoclinic lines observed with light polarized in the x direction. Int. Journ. o[ Fracture, 13 (1977) 443-454

Analysis o[ fast [racture and crack arrest by finite differences

447

Figure 3. The analytic principal stress difference field compared with a fringe pattern observed in a circular polariscope.

mesh size as small as H[40. These results are then compared with solutions corresponding to coarser grids. It was found that up to and including a mesh size of HI15, except for the crack end shape, in any other respect, no significant differences between the solutions could be detected. In particular, the stress intensity factors differ by less than 1%; the stress field maps drawings mentioned above are the same and, excluding the crack tip region, the contours of the strained specimen as plotted by the computer (Figs. lb and lc) overlap each other. Consequently, in all the dynamic analyses described below, a mesh size of h = HI15 was used. The problem has been formulated so that all physical quantities are given in terms of the specimen height H, the dilatational wave velocity C1 and the material density p. A value of 1 is substituted in the scheme for these last mentioned constants. For the problem to be completely defined (see Eqns. (1), (2) and (3)), the values of L, a0, 8 and C2 should be provided. Given the specimen geometry, L and a0 are defined in terms of H. The deflection 8, which, due to the linearity of the equations involved, may be arbitrarily chosen, is taken to be equal to H. The value of C2, however, depends on the Poisson ratio ~, of the material considered and is defined in terms of C~ according to the relation C 2 = CI(1- 2v)v2/(2 - 2~')1/2. This means that for any given geometry, different materials are distinguished in the scheme by their corresponding different Poisson ratios only. With the objective of comparing the numerical solution with experimental studies which were conducted on steel [8] and epoxy resin [7], the scheme was run, in turn, with Poisson ratio values of 0.3 and 0.392 respectively. Corresponding to the particular dimensions of the specimens used in the experiments, which happen to be similar in both studies (H[ao[L=63.5[67.81305 in [8] and 63.5[691310 in [7]), the specimen geometry used in the dynamic analysis, the same for both materials, was accordingly HlaolL = 1[1.0714.8. Shorter specimens (see Figs. 2 and 3) were used in the Int. Journ. of Fracture, 13 (1977) 443-454

448

M. Shmuely

above discussed static analysis with the benefit of saving computing time without losing accuracy since the far end of the specimen would, anyhow, not affect the static results.

3. Results for the fast crack propagation The simulation of a fracture process in the DCB runs as follows. Given the Poisson ratio of the material considered and assuming the load to be applied as illustrated in Fig. lb, we first solve the static problem. This is accompanied by solving (1) (with ot > 0), subjected to conditions (2) and (3). The solution is assumed to approximate the static state when, at any grid point, the displacement variation between two successive iterations does not exceed 10-4 8. At this stage the tryy stress level at the mesh point lying next to and ahead of the crack tip (on y = 0), denoted hereafter as crq, is recorded; the parameter ot is reduced to zero and the crack is extended by one mesh interval h. As a result of the two last mentioned changes, a dynamic state of stresses is built up in the specimen. In particular, the Cryystress level at the mesh point lying next to the tip of the just formed longer crack begins to rise. When this stress level reaches a prescribed critical value trc, the crack is further extended by one grid interval. The crack now continues to grow by one grid interval at a time as long as the tryy stress level ahead of the propagating crack tip reaches the value ~rc. The crack is assumed to be arrested when a relatively long time (t > 20 h/C2) passes after the last extension, with the tryy stress level ahead of its tip remaining below the critical value trc. During the propagation stage the time elapsing between two sequential extensions is recorded and the corresponding momentary crack velocity is evaluated. For each of the two materials investigated, the above procedure was conducted with various values of trc, each time starting from exactly the same static solution. Examples of typical cases are given in Figs. 4 and 5. Both figures illustrate that while crack propagation does not occur at a precisely constant velocity, the crack lengthtime plots, excluding the initial points, are close to being linear so that with each trq/trc ratio there can be associated a definite average velocity V. The fluctuations about this velocity may be attributed to stress waves reflected from the boundaries. In fact, reasonable correlation was obtained between the time intervals elapsing between successive points indicating velocity changes and time intervals required by the dilatational and distortional waves to span the space between the crack tip and the specimen boundaries. Both figures indicate that the crack length at arrest a= and the crack speed V increase with trq/trc. This was verified with other values of o-JtT~ as

STEEL u=0.300/ ~ . I o , , .f_~~' a- -r -r e s t ~~ 10 100 200

,_ 40

~30

o

u

0

0

Figure 4. Assorted crack length time data for two trq/trc cases in steel.

Int. Journ. of Fracture, 13 (1977) 443-454

300

Analysis of fast fracture and crack arrest by finite differences

449

,-- 401

d I

30 -

EPOXY u=0,392 -

"I-

I,-

0

i

i

I

100

200

300

t c+

-h--

Figure 5. Assorted crack length time data for two trJtrc cases in epoxy resin.

well. Notice also the sensitivity of the numerical solution to changes in the Poisson ratio which is the only parameter adjusted in the scheme when passing from one material to the other. In the course of crack propagation, at each time step, both the momentary kinetic and strain energy possessed by the specimen are calculated and their sum derived. Examples of these calculations are given in Figs. 6 and 7 which are related respectively to the cases appearing in Figs. 4 and 5. It was encouraging to find that although no energy considerations have been introduced in the analysis, it does predict a continuous decrease in the total energy as expected. Although these results cannot by themselves approve or disapprove the validity of the physical assumptions proposed here concerning the characterization of crack

1.0 =,=0.300

~3

~ '" ""

(z: 0.8

w

]--

0.6

STEEL ~ =1.5

',

~

"

mzUJ0.4--t

strain, energy

0.2

0/ 0

inetic energy

I

10

I

I

20 30 CRACK GROWTH (a - ao Yh

I

40

Figure 6. Variation of kinetic, strain and total energies in course of crack propagation for steel.

Int. Journ. of Fracture, 13 (1977) 443-454

M. Shmuely

450

v=0.392 --

EP..~_OXY

~=2.2 G~=1.6

tY hi Z tO

. 0.6

_<

,~,,,,~, / ~ t o,total t a t ener¢gy

F-

_z

0.4 ¢,r w z bA

Denergy 112

0

10

20 CRACK GROWTH

30 (a- ao)/h

40

Figure 7. Variation of kinetic, strain and total energies in course of crack propagation for epoxy resin.

growth and arrest, they do, however, indicate that at least with respect to the energetics of the system, no contradictions are involved in applying these assumptions. They also corroborate the validity of the finite difference approximation methods employed, so far established by the static solution only. Two additional important aspects of the inertial effects are emphasized by the figures. First, we notice that a considerable part of the available energy is converted into kinetic energy, most of which remains in the specimen at arrest, and should thus be taken into account in any attempt to derive an apparent dynamic stress intensity factor Ka based on the measurable initial and arrest stress intensity factors (denoted hereafter by Kq and Ka respectively) [8]. Next, the total energy curves, by their rate. of decrease, illustrate the complicated pattern of change of the energy release rate G, available for the fast propagating crack, in contrast to the constant G which would have been obtained, under quasi-static conditions, in view of the assumed constancy of the prevailing stress intensity factor K, to which the constant trc is proportional. Since, practically speaking, the velocity dependent fracture energy R = R(V) will never exactly follow the pattern of the varying G as dictated by a constant K, the actual stress intensity factor currently adjusts itself so as to provide the momentary R required. We may, however, as a first approximation, consider an apparent constant dynamic stress intensity factor Ka, which may attain different values for different trJtrc cases, as the average of the varying K-s occurring at each case, and check the consequences of this approach. The energy balance would then also be kept in the average, namely: aa

aa

fG(a)da=fR[V(a)]da. a0

ao

With this in mind, we turn to Figs. 8 and 9 where the numerical results derived Int. Journ. o[ Fracture, 13 (1977) 443-454

Analysis o[ fast [racture and crack arrest by finite differences 0.3

I

G

01

_

451

0.15 G

~- 0.2 U o._1

0.10

W

0.1

""

095

I1: U

Oq -

/

4.5

I

I

I

, ~ numerical - - , r - solution

_

4.0

I

~

-

/

experimental results

3.5

(a)

STEEL =,=0.300

\ Kq/Ktc

3.0

>. I,,Z tO I-

2.5

z

KqI ~"q'qKo

~2.0 O: U.I

.J 1.5 U.l

,.0

(b) J

1.oL~" 1

I

2

I

3

RELATIVE ARREST LENGTH

t, a=/ao

Figure 8. Relation between crack growth, crack velocities and stress intensity factors in steel.

here for steel and epoxy resin are compared with experimental measurements reported in [8] and [7] respectively. In Fig. 8 the experimental plots of V/C2 and were enlarged from Fig. 4c of [8]. From the Kq/~/K~/K= data, together with the information provided in Fig. 8 of [8], the KqlK~c values were derived. The variables V, Kq and K= are respectively the measured average velocity and the calculated initial and arrest stress intensity factors, and Kit is the mode I fracture toughness of the material considered. By the numerically calculated crack velocity in Fig. 8a, we mean the constant value least square approximation of the velocities encountered in any particular trJ~c run.

KJ'k'/~Ka

Int. Journ. of Fracture, 13 (1977) 443--454

M. Shmuely

452 0.3

0.2

(a)

1

O.05

lexp~mentol ~) numerical EPOXY u=0.392

0.1

7 nO

0 3.0 _

I

K~xo

_-

u')

~2.5 0

_

I~IK~

~ 2.0 t5

(b)

1.0 1

2

3

4

RELATIVEARREST LENGTH aolao Figure 9. Relation between crack growth, crack velocities and stress intensity factors in epoxy resin.

Both the numerically evaluated crack speed V and the crack length at arrest have been found to be single valued increasing functions of trJtrc as illustrated in Figs. 8 and 9. This implies a unique relation between V and a,, which is the only measurable relation which can be compared in steel. With the good matching (Fig. 8a) obtained in hand, given the experimentally measured Kq and the analytically assumed o'q/trc as related to either aa (Fig. 8b), or alternatively to V, the apparent dynamic stress intensity factor can be determined according to the relation

(4) or

K~ = 5 Ka 0% As argued in [8], Kd may be alternatively approximated by the relation Int. Journ. of Fracture, 13 (1977) 443-454

(5) K d =

~v/KqKa,

Analysis o[ fast fracture and crack arrest by finite differences

453

which is a close approximation when all the kinetic energy is recovered at arrest. Following relation (5), the KJX/KqKa curve is compared with the trq/trc curve to find that the first variable deviates from the other by an almost constant factor (notice that a logarithmic scale is used) of about 1.06, which may be attributed to neglecting the unrecovered kinetic energy in the Ka approximation. Part of the experimental data obtained and reported in the epoxy resin study [7] was different from that corresponding to steel. Using the shadow spot method of Manogg in combination with a Cranz-Schardin high speed camera, it was possible to trace and determine the current dynamic stress intensity factors Kd and the current crack speeds V. In all cases tested, a steady state velocity and a corresponding almost constant Ka were observed during most of the propagation stage with both variables gradually decreasing towards the arrest. In Fig. 9a the steady state velocities rather than the average velocities are depicted. While the lower curve of the numerical results in this figure was obtained in a similar manner as that previously described for steel, in the upper curve only the first 2/3 of the velocities encountered in a given run were considered in the least square approximation, with the objective of considering that part of the propagation stage to which the experimental results are related. Note that here again a good matching is obtained. In the KJKd curve, rather than consider an approximated Kd as in steel, we substitute for Kd the weighted average of the Kd-S observed in the experiments. Again, we find a good agreement between the KJKd values and the corresponding trq/trc values with deviations less than 7%. Comparing the relative levels of the Kq/K~ and trJ~rc curves in both figures and by referring to relation (4), we find that in terms of K~c, the dynamic stress intensity factors in steel are in general higher than those encountered in epoxy resin when the same relative velocities or crack growth are considered. Notice also that for very small velocities the KJKt~ curve tends to be (Fig. 8) or actually is (Fig. 9) below the trq/trc curve which means that for these velocities K d is smaller than K~. With the detailed information on Kd as provided in [7], it would be only natural to consider an improved simulation where the tr~ in the scheme constantly adjusts itself to values implied by the instantaneous velocities encountered. Also, the effect of changes in geometry and material is of great interest. These matters will be studied in the near future.

REFERENCES [1] M. Shmuely and Z.S. Alterman, Journal of Applied Mechanics, 40, 4, Transactions of ASME, 95, Series E (December 1973) 902-908. [2] M. Shmuely and D. Peretz, International Journal of Solids and Structures, 12, 1 (1976) 67-79. [3] K.B. Broberg, "The Propagation of a Brittle Crack", Arkiv [6r Fysik, Vol. 18, Stockholm, Sweden (1%0) 15%192. [4] J.R. Rice, "Mathematical Analysis in Mechanics of Fracture", in Fracture, Vol. 2, ed. Liebowitz, H., Academic Press (1%8) 191-213. [5] L.B. Freund, Journal of the Mechanics and Physics o[ Solids, 20 (1972) 12%152. [6] M.F. Kanninen, "An Analysis of Dynamic Crack Propagation and Arrest for a Material Having a Crack Speed Dependent Fracture Toughness", in Prospects of Fracture Mechanics, ed. Sih, G.C. et al., Noordhoff, Leyden, The Netherlands (1974). [7] J.F. Kalthoff, S. Winkler and J. Beinert, International Journal of Fracture, 12 (1976) 317-19 and private communication. [8] G.T. Hahn, R.G. Hoagland, M.F. Kanninen, A.R. Rosenfield and R. Sejnoha, "Fast Fracture Resistance and Crack Arrest in Structural Steels", SSC-242 Progress Report on Project SR-201, "Fracture Arrest Study", Dept. of the Navy (1973) 1-24. Int. Journ. of Fracture, 13 (1977) 443-454

454

M. Shmuely

RI~SUMI~ Un sch6ma par diff6rence finie propos6 r6cemment et relatif principalement h la solution statique d'une 6prouvette double cantilever sous conditions d'6tat plan de d6formation, est 6tendu h I'obtention d'informations compl6mentaires sur les conditions de propagation instable et d'arr6t d'une fissure. L'analyse pr6dit des Iongueurs de fissure h I'arr~t et des vitesses de fissuration, ainsi que les facteurs d'intensit6 de contrainte dynamique qui sont tr~s proches de ceux que r o n trouve dans la litt6rature. Ce bon accord avec les observations exp~rimentales a 6t6 obtenu pour de l'acier et pour une r~sine 6poxy, deux mat6riaux qui, bien que de classes diff6rentes, ne se distinguent n6anmoins dans le cadre de ranalyse que par leur modules de poissOn.

Int. Journ. of Fracture, 13 (1977) 443--454