Dynamic Fracture
5.6 Examples for Propagating Cracks
1− t
t∗
, (5.100)
the correction term associated to the power of order5/2can be calculated using the DENG’s solution and by setting for exampleb51=b52= 0. The particular values KImax = 10 Kgs2/mm1/2 andt∗ = 50μswere used. The body under study is the same as in the convergence analysis in Section 4.5.3, that is, the unit squareΩ = [0, 1]×[0,1]of side length1 mcontaining a crack going from(0,0.5)to(0.5,0.5). The material parameters areE= 200 Pa,ν= 0.3andρ= 5000 Kg/m3, which together with a plane strain model give the wave velocitiesCd= 7338 m/s,Cs= 3922 m/sand CR= 3638 m/s.
The problem is solved using the Generalized-αmethod withρ∞= 1and1280time steps. Furthermore, three meshes of25×25,45×45and75×75bilinear quadrilateral elements are considered. The SIFKIis extracted using a fixed influence radius of 0.2 m.
The solutions are shown in Figure 5.21 and compared with the exact solution (5.100). From this figure, it is clear that the oscillations of the numerical solutions are an artefact of the spatial discretization and they will vanish as the mesh size tends to zero. Another parameter which influences these oscillations ist∗; small val-ues oft∗ produce a fast variation ofKI, which is difficult to capture with a coarse mesh, as shown in Figure 5.21. On the other hand, bigger values oft∗generate a slower variation ofKIwhich can be described even with a coarse mesh. This is illustrated in Figure 5.22 where the valuet∗= 200μsis chosen. In this case even the 25×25mesh can capture the behaviour ofKIwith good accuracy.
5.6 Examples for Propagating Cracks
Turning now the attention to a propagating crack. The main problem in modelling this phenomenon is that the enrichment changes as the crack propagates. Thus, the problem can be formulated as follows: Given the displacement, velocity and acceleration at timetn(un,u˙nandu¨n), associated with the old enrichment, together with the crack tip position, crack velocity, and crack propagation angle (an,vnand θn) find:
• un,u˙nandu¨n, at timetnbut associated with the new enrichments.
• un+1,u˙n+1andu¨n+1.
• an+1,vn+1andθn+1.
0 0.1 0.2 0.3 0.4 t/t∗ 0.6 0.7 0.8 0.9 1.0
−0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2
Mesh25×25 Mesh45×45 Mesh75×75 u
Exact (5.100) KI/KImax
Figure 5.21:KI/KmaxI using the Generalized-αmethod (ρ∞= 1) and three meshes.
t∗= 50μsandΔt=t∗/1280.
0 0.1 0.2 0.3 0.4 t/t∗ 0.6 0.7 0.8 0.9 1.0
−0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2
Mesh25×25 Mesh45×45 Mesh75×75 Exact (5.100) KI/KImax
Figure 5.22:KI/KmaxI using the Generalized-αmethod (ρ∞= 1) and three meshes.
t∗= 200μsandΔt=t∗/1280.
In order to solve this problem,an+1is calculated explicitly by adding a new crack segment of lengthΔa = vnΔtin the direction of the crack propagation angleθn, thus the new crack geometry is known. Next, two strategies are considered for the determinationun,u˙nandu¨n:
1. As proposed in [Réthoré et al., 2005b], the old singular enrichments are kept (of course associated with its respectiveoldtip) and new singular enrichments are added (associated with the new crack tip). The discontinuous enrichments are selected in the usual way. The idea is that all the old enrichments are present in the new set of enrichments, and hence the natural choice is to ini-tialize the new enrichments to zero, that is,
un=
uTn|0 0...0 0T
, u˙n=
u˙Tn|0 0...0 0T
, u¨n=
u¨Tn|0 0...0 0T
.
(5.101)
It has been proved in [Réthoré et al., 2005b] that this scheme introduces no numerical energy due to the evolving enrichments. From now on, this method will be referred to as ECS (Energy-Conserving Scheme). See also [Elguedj et al., 2009] for some examples of using this technique in conjunc-tion with an explicit time integraconjunc-tion algorithm and [Grégoire et al., 2007] for a comparison with experiments, considering the NEWMARK’s implicit mean ac-celeration scheme as time integrator.
2. A different strategy is proposed in this thesis in the framework of the XFEM, which consists in determining the new enrichments in the usual manner (with-out keeping the old singular enrichments) and then mappingun(and alsou˙n
andu¨n) to the new enrichment configuration by minimizing the error e=
Ω
(un(X)−un(X))·(un(X)−un(X))dΩ, (5.102) which leads to
Anun=Anun, (5.103)
with An=
Ω
(Φn)TΦndΩ and An=
Ω
(Φn)TΦndΩ, (5.104) whereΦnandΦnare the matrices containing the shape functions at timestn
for the old and new crack configurations, respectively, that is, Φn= [φ1I3x3 φ2I3x3 φ3I3x3...φnbfI3x3],
Φn=
φ1I3x3 φ2I3x3 φ3I3x3...φnbfI3x3 ,
(5.105)
beingnbfandnbfthe number of XFEM basis functions in each configuration.
The mapping ofu˙n andu¨n follows straightforwardly from (5.103) using the same matricesΦn andΦn. The method will be referred to as MS (Mapping Scheme).
Once the displacement, velocity and acceleration are calculated for the new enrich-ments at timetn, they can be calculated at timetn+1by using one of the time integra-tion schemes described in secintegra-tions 5.3 and 5.4. At this point, it is worth menintegra-tioning that since the acceleration is prescribed, the crack forces (Fcrack) necessary to close the crack at timetnfor the new enrichments must be calculated, for example, using the equation of motion
Mnu¨n+Dnu˙n+Knun=r(tn) +Fcrack(tn). (5.106) Of course, if the NEWMARKmethod is used, then the crack forces at timetnplay no role and there is no necessity of using 5.106. On the other hand, if the dG method is considered, then the time dependency ofFcrackover the time interval[tn, tn+1]is needed. In this case,Fcrackis approximated linearly over the time interval, that is,
Fcrack(t) =Fcrack(tn)(1−(t−tn)/(tn+1−tn)). (5.107) Onceun+1,u˙n+1andu¨n+1are known, the stress intensity factors can be extracted using the interaction integral (Section 3.6). Knowing the SIFs, the crack initiation, crack growth and crack arrest criteria reviewed in Section 3.7 can be used with one of the mixed-mode dynamic crack propagation criteria discussed in sections 3.8, 3.9 and 3.10 in order to evaluate whether the crack propagates or not, and if yes, to calculate the crack propagation angleθn+1and crack tip velocityvn+1. For example, using the maximum circumferential stress criterion, the equivalent stress intensity factorKeqcan be used to determinevn+1as follows
vn+1=
⎧⎪
⎪⎪
⎪⎪
⎪⎪
⎨
⎪⎪
⎪⎪
⎪⎪
⎪⎩
0 if Keq< KId and vn= 0, v0(1−KIA/Keq)1/m if Keq≥KId and vn= 0, 0 if Keq< KIA and vn>0, v0(1−KIA/Keq)1/m if Keq≥KIA and vn>0.
(5.108)
Here, theK−vrelationship given in (3.98) was considered, wherev0andmare empirical parameters,KIAis the crack arrest toughness andKIdis the crack initiation toughness.
5.6.1 Semi-infinite Crack in an Infinite Plane Subjected to a Stress Pulse
This problem has been solved analytically by FREUND [Freund, 1990] and used to test numerical methods by several authors, see [Nakamura et al., 1985, Belytschko et al., 2003, Réthoré et al., 2005b]. The material parameters of the plane strain model areE = 210 GPa,ρ= 8000 Kg/m3andν= 0.3. The geometry of the model is defined byH= 2 m,L= 10 manda= 5 m, see Figure 5.23.
For a stationary crack, the analytical solution is given by KI0(t) =
⎧⎪
⎨
⎪⎩
0 if t≤td, 2σ0
1−ν
Cd(t−td)(1−2ν)/π if t > td,
(5.109)
L
a 2H
σ0
Figure 5.23: Geometry to simulate a semi-infinite crack in a infinite plane subjected to a stress pulseσ0.
wheretdis the time of arrival of the dilatational waves to the crack tip, this is,td = H/Cd. This solution is valid until the reflected waves from the boundary reach the crack tip, at timeT = 3td. The problem was solved using a mesh of80×41bilinear quadrilateral elements, a total calculation time of3tdand 200 time steps. The SIFKI
was calculated using an influence radius of0.5 m. The results using the Generalized-αmethod withρ∞ = 1are shown in Figure 5.24, here, three different cases are considered: (a) only the nodes whose support contains the crack tip are enriched by the near tip functions (amax= 0), (b) all the nodes within a radius ofamax= 0.75 mare enriched by the near tip functions and (c) the row-sum lumped mass matrix is used for the elements without enrichment, together withamax= 0.75 m. All the solutions show a good agreement with the analytical solution, but lumping the mass matrix eliminates the unwanted oscillations before the arrival of the dilatational waves at the crack tip. Very similar results are obtained using the dG method withp= 1as shown in Figure 5.25. Another comparison is carried out using only 20 time steps, amax = 0and no lumping for the cG and dG methods considered in this chapter.
Figure 5.26 shows these results for the cG method and polynomials of degree 1, 2 and 3; the good quality of the solution using high order polynomials and only 20 time steps is evident. Similar results are obtained for the dG method as shown in Figure 5.27, but with better accuracy in comparison with the cG method, mainly in the casesp= 1andp= 2.
Now, the case of a propagating crack that starts to propagate at timet=tdwith a prescribed velocityv= 1500 m/sis analysed. The analytical stress intensity factor KI(t)is given byKI(t) = kI(v)KI0(t)(see Section 3.4.3), where for this example kI(v) = 0.652. For the calculationsamax was set to zero, no lumping was consid-ered and only 20 time steps were used. Figure 5.28 compares the analytical so-lution forv = 1500 m/swith the results of using the ECS and MS together with the Generalized-αmethod as time integration scheme. The analytical solution for v= 0is also plotted as reference. Even when the results are very similar for both schemes, the quality of the solution decreases in comparison with the results ob-tained forv= 0. Slightly better results are obtained using the dG method as shown in Figure 5.29, where also the mean value of the jumps for the MS are plotted. Again the results using the ECS and MS are very similar. A final investigation is carried out using a finer mesh around the crack propagation zone (Figure 5.30), which consists
0 0.5 1.0 t/td 2.0 2.5 3.0 KI/σ0√
H
−0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
amax= 0
Analytical solution
amax= 0.75 m amax= 0.75 mand lumpedM
Figure 5.24: Results ofKI/σ0√
Husing the Generalized-αmethod withρ∞= 1and 200 time steps.
0 0.5 1.0 t/td 2.0 2.5 3.0
KI/σ0√ H
−0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
amax= 0
Analytical solution
amax= 0.75 m amax= 0.75 mand lumpedM
Figure 5.25: Results ofKI/σ0√
Husing the dG method withp = 1and 200 time steps.
in 5054 elements. The elements in the centre have a size of2.67×2.85cm while the size of the elements near to the boundary is of21×23cm. Figure 5.31 shows the results using the mapping scheme and the dG method (p= 1), now with 40 time steps. Clearly, the quality of the solution is improved.