# Study on the Fracture Process Zone near the Mode I Dynamic Crack Tip

Authors: Chenmeng Ji; Chengzhi Qi
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Evaluation of the shape and size of the fracture process zone near the mode I dynamic crack tip is still a problem unsolved completely at present. The research on the relationship between the fracture process zone and crack velocity near the mode I dynamic crack tip is quite limited, and some researchers have also developed experimental methods or numerical methods. In this research, based on the theory of elastodynamics and the complex stress function method, an approximate method for solving the mode I dynamic crack problem was proposed. The fracture process zone near the mode I dynamic crack tip was analyzed. The results showed that the areas of the fracture process zone determined based on the approximate method are nearly the same as the results obtained based on the well-known stress fields. The approximate method could provide a good reference for determining the fracture process zone near the mode I dynamic crack tip since no analytic methods had been found for evaluating the fracture process zone near the dynamic crack tip to the authors’ knowledge.


---

## Body

## 1. Introduction

The study on the fracture process zone has been mainly focused on the fracture process zone near the static crack tip, while the analysis of the fracture process zone near the dynamic crack tip is very complicated due to the inertia and the dynamic effect of materials. For the fracture process zone near the static crack tip, some analytical models, such as the cohesion model [1], the virtual crack model [2], and two-parameter model [3], treat the fracture process zone as linear in form. The cell model proposed by Broberg [4–7] well explained the existence of the fracture process zone at the crack tip and provided a very visual explanation for future researchers. At present, the research studies on the dynamic crack tip fracture process zone are very limited, and some researchers have also developed experimental methods or numerical methods. For the determination of the fracture process zone near the dynamic crack tip, although Broberg [8] and Fineberg [9] concluded that the fracture process zone at the dynamic crack tip increases with the increase of the crack velocity, the authors have searched relevant literature and found no relevant analytical methods yet. The determination of the profile of the fracture process zone at the dynamic crack tip is still an unsolved problem. The integral transformations with dual integral equation method and Wiener–Hopf method are difficult to allow for obtaining analytical solutions due to the complexity of calculations. The complex function method is much simpler than the integral transformation methods because it does not need to carry out the complex inversion operations.In this study, based on the theory of elastodynamics and the complex stress function method, an approximate method for solving the dynamic crack problem was proposed. The fracture process zones at the mode I dynamic crack tip were analyzed. The contours and areas of the fracture process zones were determined by using the maximum tensile stress criterion for different crack velocities. Finally, the results based on the approximate method proposed in this paper were compared with the results calculated based on the well-known stress field at the crack tip to verify the feasibility and correctness of the approximate method.


## 2. Method

The crack tip is assumed to move along the direction ofx-axis at a constant speed V. The far-field tension σ is applied, see Figure 1. The fixed coordinate systems (x1, y, t) and the moving coordinate system (x, y, t) have the following relationships:(1)x=x1−Vt,y=y.Figure 1
Model of mode I crack.The coordinate transformation gives(2)∂x∂t=−V∂∂x.For plane elastic problem, the equation of motion takes the following form:(3)∂σx∂x+∂τxy∂y=ρ∂2u∂t2,∂σy∂y+∂τxy∂x=ρ∂2v∂t2,where ρ is mass density and u and v are displacements along x and y axes.The relationships between stresses and strains give [10](4)σx=λθ+2μ∂u∂x,σy=λθ+2μ∂u∂y,τxy=μ∂u∂y+∂v∂x,θ=∂u∂x+∂v∂y,where μ and λ are Lame’s constants and θ is the first invariant of strain tensor, that is,(5)2μ∂u∂x=σx−λ2λ+μσx+σy,2μ∂v∂y=σy−λ2λ+μσx+σy,θ=12λ+μσx+σy.The equation of compatibility takes the form as follows [10]:(6)∂2∂x2+∂2∂y2σx+σy−2λ+μλ+2μρ∂2∂t2∂u∂x+∂v∂y=0.Substituting equation (5) to equation (6), the following equation can be obtained:(7)∂2∂x2+∂2∂y2−ρλ+2μ∂2∂t2σx+σy=0.Taking partial derivatives with respect tox and y for both sides of equation (3), the following equations can be obtained:(8a)∂2σx∂x∂y+∂2τxy∂y2=ρ∂3u∂t2∂y,∂2σy∂x∂y+∂2τxy∂x2=ρ∂2v∂t2∂x,(8b)∂2σx∂x2+∂2τxy∂x∂y=ρ∂3u∂t2∂x,∂2σy∂y2+∂2τxy∂x∂y=ρ∂3v∂t2∂y.By adding and subtracting both sides of the above two groups of equations, the following can be obtained:(9)∂2σx∂x2−∂2σy∂y2=ρ2μ∂2∂t2σx−σy.Rearranging equation (9) yields(10)∂2∂x2−ρ2μ∂2∂t2σx=∂2∂y2−ρ2μ∂2∂t2σy.Obviously, introducing the stress functionU yields the following equations [10]:(11)σx=∂2∂y2−ρ2μ∂2∂t2U,σy=∂2∂x2−ρ2μ∂2∂t2U.Adding the above two equations yields the following equation:(12)σx+σy=∂2∂y2+∂2∂x2−ρμ∂2∂t2U.Substituting equation (12) into equation (7) yields(13)∂2∂x2+∂2∂y2−ρλ+2μ∂2∂t2∂2∂x2+∂2∂y2−ρμ∂2∂t2U=0.By substituting some variables, equation (13) is transformed as follows:(14)∂2∂x2+∂2∂y2−1c12∂2∂t2∂2∂x2+∂2∂y2−1c22∂2∂t2U=0,where c1 is the longitudinal wave velocity and c2 is transverse wave velocity, which are defined as follows:(15)c1=λ+2μρ,c2=μρ.According to equations (2) and (14), it can be changed to the form as follows:(16)∂2∂x2+∂2∂y12∂2∂x2+∂2∂y22U=0,where(17)α1=1−V2c12,α2=1−V2c22,y1=α1y,y2=α2y.Since the analytical solution of equation (16) cannot be obtained accurately, and the expression of the stress function of the dynamic crack cannot be determined, the stress field near the mode I crack tip cannot be determined directly by solving equation (16). The controlling equation of the static problem is one biharmonic equation, and the stress function method for solving static fracture process zone near the crack tip is already very complicated. Equation (16) is very similar in form to the biharmonic equation. Hence, we consider using the method for solving static crack problems to solve the dynamic problem of the stress field near the crack tip.Let(18)y¯=21/α12+1/α22y=2α12α22α12+α22y.So, equation (16) can be rearranged as follows:(19)∂4∂x4+2∂4∂x2y¯2+∂4∂y¯4U+4α12α22α12+α222−1∂4∂y¯4U=0.Let(20)β=4α12α22α12+α222−1.So, equation (19) can be rearranged as follows:(21)∂2∂x2+∂2∂y¯22U+β∂4∂y¯4U=0.It can be seen that equation (21) is a fourth-order partial differential equation about x and y¯. It is very difficult to solve equation (21) without some efforts. However, the second term on the left side of the equation contains parameter β, and the values of β at different crack velocity are shown in Table 1.Table 1
The values ofβ for different velocities.
V/cα1α2α12α22β0111100.190.99450.99840.98910.9968−1.53E − 050.380.97800.99370.95650.9875−0.000250.570.94980.98590.90210.972−0.001390.760.90880.97480.82600.9503−0.004900.950.85330.96040.72810.9224−0.01386It can be seen from Table1 that β increases with the crack velocity. The maximum value of β is reached when the Rayleigh wave velocity is approached. However, the generally observed maximum crack propagation velocity in experiments is 0.34–0.45 times of the Rayleigh wave velocity, and the value of β is in the range of −0.00025 to −0.00139. The value of β is very small, and the second term of equation (21) can be neglected. Hence, the following equation can be used to explore the approximate solution of equation (16):(22)∂2∂x2+∂2∂y¯22U=0.At this time, equation (15) is the biharmonic equation in the coordinate system (xy¯) and has the same form as that for the static crack. Then, we can find the approximate solution for the dynamic crack by using the stress function method.For mode I crack, Westergaard [11] assumed the stress function taking the following form:(23)Ux,y=Re∫dz∫Zzdz+yIm∫Zzdz.The complex stress function is [12, 13](24)Zz=K12πz11+34za−532za2+21384za3−….Then, the stress fields near the mode I crack tip can be calculated:(25)σx=ReZ−y¯ImZ′,σy=ReZ+y¯ImZ′,τxy=−y¯ReZ′.The size of the FPZ can be calculated theoretically based on the appropriate criterion and the crack tip stress field. The maximum tensile stress criterion for brittle materials can be expressed as follows:(26)σ1=σx+σy2+σx−σy22+τxy2=σu,where σu is the tensile strength and is assumed to be independent of the crack velocity.When more than one term in the expansion of the stress function is taken, the expressions for the stress field based on the stress function method at the crack tip are shown in Appendix.


## 3. Results

By substituting equation (25) into equation (26), the implicit equation of the FPZ can be easily obtained. As an example, we take the following material parameters: the longitudinal wave velocity c1 = 5370 m/s, the shear wave velocity c2 = 3180 m/s, and the Rayleigh wave velocity c = 2900 m/s. The contours of the FPZ constructed by taking different terms of the stress function of equation (24) for different crack velocities are shown in Figure 2.Figure 2
The FPZs constructed by using the maximum tensile stress criterion with different terms of stress function: (a) only first term, (b) first two terms, (c) first three terms, and (d) first four terms.
(a)(b)(c)(d)From Figure2, it can be seen that the fracture process zone increases with the increasing crack velocity. The higher the velocity is, the more significant the change of the fracture process zone is. The area of the fracture process zone at 0.38c is 1.026 times that at 0.19c, and the area of the fracture process zone at 0.76c is 1.34 times that at 0.95c, or 1.657 times that at 0.19c. Taking the velocity of 0.38c as an example, the calculated areas of the fracture process zone are shown in Table 2. The area of the fracture process zone calculated by the first two terms of the stress function is the largest. The more the stress function terms are taken, the smaller the change rate of the area of the fracture process zone at the crack tip is, which also indicates that the calculation results become more and more stable.Table 2
Areas of the FPZ for different terms in the expansion of the stress function.
Stress functionArea of the fracture process zone (cm2)Change rate (%)First term4.23987—First two terms7.8160484.35First three terms7.314736.41First four terms7.326980.167First five terms7.323140.052


## 4. Discussion

The approximate method is proposed in the second part by ignoring the influence of the less-significant parameters. The well-known stress field was used to calculate the fracture process zone at the crack tip, and the results by using the approximate method were compared to the results by using the well-known stress fields to verify the correctness of the approximate method. The expressions for the mode I stress fields are given as follows [14]:(27)σx=∑Kn0BIV2πnn+121+2α12−α22r1n/2−1cosn−1θ12−2hnr2n/2−1cosn−1θ22,σy=∑Kn0BIV2πnn+12−1+α22r1n/2−1cosn−1θ12+2hnr2n/2−1cosn−1θ22,τxy=∑Kn0BIV2πnn+12−2α1r1n/2−1sinn−1θ12+1+α22α2hnr2n/2−1sinn−1θ22,where(28)BIV=1+α224α1α2−1+α222,hn=2α1α21+α22,odd,1+α222,even,x¯=x1−a˙t=r¯1cosθ1=r¯2cosθ2,y¯=y=r1asinθ1=r2asinθ2,r1=x¯2+α12y¯2,r2=x¯2+α22y¯2,θ1=arctanα1y¯x¯,θ2=arctanα2y¯x¯.If equation (27) is substituted into equation (26), we can obtain the equation for the fracture process zone at the crack tip r = f (θ). The higher-order terms of the stress field have little effect on the stress field near the crack tip of mode I, but the second term has a significant effect on the fracture process zone at the crack tip. Therefore, the first term and the first two terms in the stress field expansion expression are considered, respectively, to calculate the fracture process zone at the crack tip of mode I. K20 is a parameter similar to the stress intensity factor, which is difficult to determine. Because the dependence of T-stress on the velocity for dynamic crack is not clear, the ratio of T-stress to the far-field stress is 0.526 referred to the work by Seed and Nowell [15].When the first term and the first two terms of the stress function expansion are taken respectively, the comparisons of the areas of the fracture process zone at the mode I crack tip calculated based on the approximate method and the well-known stress field expression are shown in Table3. It can be seen that the results of the mode I crack tip fracture process zone calculated by the approximate method are highly consistent with those of the crack tip fracture process zone calculated based on the known stress field, and both increase with the increase of the crack velocity. When the velocity is relatively high, the area of the fracture process zone at the mode I crack tip calculated by the approximate method is slightly smaller than that calculated based on the well-known stress field, which is believed to be caused by neglecting the influence of small parameters β. The comparison of the calculation results of the two methods in Table 3 verifies the correctness of the approximate method in this paper for evaluating the fracture process zone at the tip of mode I dynamic crack, which can provide a reference for the study of the mode I dynamic crack problem.Table 3
Comparison of the FPZ areas by using the approximate method and the method based on the well-known stress (cm2).
Crack velocity (m/s)1003005007009001100First termWell-known stress method4.09474.16454.26694.42684.69405.0909Approximate method3.97804.00314.05784.09304.14004.2102First two termsWell-known stress method6.22646.37386.60777.00347.60298.5297Approximate method7.46947.47747.50427.57127.66787.7990


## 5. Conclusions

An approximate method was proposed by transforming the governing equation and disregarding the less-significant parameter. The areas and sizes of the fracture process zone of the mode I dynamic crack were evaluated by the approximate method. The correctness of the method was verified by comparing the results calculated by the approximate method and the well-known stress fields. The conclusions can be drawn as follows:(1)
The fracture process zone near the mode I dynamic crack tip distributed symmetrically about the crack surface. The higher the crack propagation velocity, the larger the size of the FPZ for the mode I crack. When the velocity approached the Rayleigh wave velocity, the fracture process zone changed more significantly.(2)
The fracture process zone can be calculated by using this approximate method with different terms of the stress function. The sizes and areas of the fracture process zone calculated by the proposed method were very close to those calculated based on the well-known stress field. Therefore, the proposed method is suitable for the mode I dynamic crack.(3)
Although the well-known stress field is widely accepted, it is just an asymptotic solution for the mode I dynamic crack tip, i.e.,n = 1. When n > 1, the parameter K20 cannot be determined except by experiment. This means that it is difficult to use the well-known stress field expression to evaluate the FPZ for the dynamic crack tip for n > 1. However, the approximate method can be used to obtain the results with different accuracies by taking different terms from the expansion of the stress function. It is also noteworthy that the method is very easy to implement for determining the shape and size of the FPZ of the mode I dynamic crack.



---

*Source: 5515004-2021-05-03.xml*