XFEM for Crack Modelling
4.4 Numerical Aspects
4.4.2 Node Selection for Enrichment in Order to Avoid Linear DependencyDependency
The node selection for the enrichments using the HEAVISIDEfunction previously de-scribed (discontinuous enrichment, intersection enrichment and branching enrich-ment) needs a further step in order to avoid linear dependency. To clarify this, con-sider a node, nodei, having the discontinuous enrichment (HEAVISIDE), see Figure 4.23. Suppose, without loss of generality, that the functionH(X)takes the value 1 on the side of the crack containing the node under study. Then, the functionsφi(X) andφi(X)H(X)differ only on the areaA2, represented in the Figure 4.23 b). Now, if this area is too small compared with the area of the nodal supportAw, there is only a
little difference between the two functions and the nodeiis enriched practically with the same functionφi(X), which would lead to linear dependency. The same linear dependency problem appears ifA1is too small. To avoid this situation, the size of the areas is used to define the enrichment. Formally, if occurs that
min (A1, A2) Aw
< tol, (4.75)
the node is removed from the set of nodes with HEAVISIDE enrichment. Dolbow [Dolbow, 1999] proposed a tolerance of0.01%, which has shown good results for fracture problems.
Crack
A1
A2 Support
Node under study
a) b)
Figure 4.23: Node selection strategy to avoid linear dependency: HEAVISIDE enrich-ment.
The cases of intersection and branching are treated in a similar way. Consider Figure 4.24, in this case the node under study is eliminated from the set of nodes with intersection enrichmentI(X)if
min (A1, A2, A3, A4) Aw
< tol. (4.76)
Analogously, the node in Figure 4.25 is eliminated from the set of nodes with branch-ing enrichmentJ(X)if
min (A1, A2, A3) Aw
< tol. (4.77)
The sub-areas can be easily calculated using the triangulation described in the pre-vious section.
4.4.3 Transformation for a kinked crack Considering that the near-tip functionF1=√
rsinθ/2allows to model the disconti-nuity along the crack due to the different values of the function forθ=πandθ=−π.
Crack1 Crac
k2
A1
A2 A3 A4
Support Node under study
a) b)
Figure 4.24: Strategy to avoid linear dependency: Intersection enrichment.
Crack1 Crack 1
Crac k2
A1
A2 A3
Support Node under study
a) b)
Figure 4.25: Strategy to avoid linear dependency: Branching enrichment.
Then, for kinked cracks, a mapping is needed to align the crack segments to the leading segment, in order to haveθ=±πon the crack segments influenced by the near-tip functions.
The mappings starts by defining the angleαref, which determines a reference line that divides the domain into two regions, see Figure 4.26. On the region contain-ing the leadcontain-ing segment no transformation is considered, while on the other re-gion the transformation takes actually place. The angleαref can be simply defined as a constant, in this case the usual choice isαref = π/2([Fleming et al., 1997], [Belytschko and Black, 1999]), or by the bisector of the angle formed by the vectors X2−XkinkandXtip−Xkink. In one way or another,αrefcan be considered known and being dependent at most on the crack geometry.
The idea is now to transform the points on the left of the reference line, in such a way that points over this line are not affected but points over the secondary crack segmentX2Xkinkare aligned to the leading segment. A simple way to perform this, it is to rotate the points on the mapping zone about the kink pointXkink. Thus, ifβ
e2
Figure 4.26: Mapping for a kinked crack.
is defined as the angle between the vectorv=X−Xkinkand the reference line (measured anticlockwise), andβrotis the rotation angle of the vectorvabout the kink pointXkink, the following conditions need to be satisfied:
βrot=
A simple way to satisfy these conditions is to carried out a linear rotation, as in [Fleming et al., 1997], imposing continuity over the reference line and the secondary segment, formally
With this definition, the functionβrot(β)is only of classC0, due to the discontinuity of its derivative atβ= 0andβ =−π. Now, if the fact that the not only the
map-ping but also its derivative are involved in the determination of the stiffness matrix is considered, the discontinuity in the derivative would require, from the theoretical point of view, a modification of the triangulation exposed in 4.4.1 to take the discon-tinuity into account. In order to avoid this, a mapping of classC1over the domain β∈[−π, α)∪(α, π]can be defined as
Once the rotation angle is calculated, the rotation matrixMcan be defined as
M=
These new coordinates are used to calculaterandθfor the evaluation of the near-tip functions, with
Since not only the coordinatesrandθbut also its derivatives with respect to the position,∂(r, θ)/∂X, are needed. The application of the chain rule leads to
∂(r, θ)
and the partial derivative∂X/∂Xcan be calculated from (4.82), resulting in
∂X
Equation (4.87) can be written as
∂X
∂X =M+∂βrot
∂β
∂M
∂βrot
X−Xkink ∂β
∂X T
(4.89) where∂βrot/∂βis obtained from (4.80) (or (4.79) if the reader insists in using theC0 mapping),∂M/∂βrotcan be calculated from (4.81) and∂β/∂Xis given by
∂β
∂X = 1
X−Xkink2
⎡
⎣−
X2−X2kink X1−X1kink
⎤
⎦. (4.90)
Figure 4.4.3 shows the functionF1for a crack with a kink of45◦and different mapping alternatives. The function is plotted over the domain [−1.5 0.5]×[−1 1]and the mappings considered are:
• Mapping of classC0withαref= 90◦, Figure 4.4.3(b).
• Mapping of classC1withαref= 90◦, Figure 4.4.3(c).
• Mapping of classC1withαref = 112.5◦, which corresponds to the bisector of the angle at the kink point, Figure 4.4.3(d).
Figures 4.28 and 4.29 show the functionF1and its partial derivative∂F1/∂X1for the different mappings as function ofX1, settingX2 =−1.0. In the first figure, the kink in the plot of theC0mapping at the pointX1=−0.5(on the reference line) can be recognized. This kink implies a discontinuity in the partial derivative∂F1/∂X1at the reference line, as shown in Figure 4.29 forX2=−1.0.
In the case of a crack with several kinks, the mapping described above can be applied successively in order to determine the final coordinatesrandθfor the eval-uation of the near-tip functions. Consider the arbitrary pointXin Figure 4.30. The first step is to determine how many mappings are needed. This can be carried out by calculating the angleβassociated to the point Xfor each kink point, starting from the leading segment. If for a kink pointβ∈[0,−π], it means that no mapping with respect to that kink is necessary and the mappings start from the previous kink point. This is schematized in Figure 4.30(a) where the mappings are needed from the second kink point. The successive mappings are also shown in the figure. First the pointXis mapped intoX(1)by rotating aboutXkink(1) and the third segment is aligned with the second segment. The result of the first mapping is shown in Fig-ure 4.30(b). Next, a new mapping is applied, rotatingX(1)aboutXkink(2) and aligning the three segments. Thus, the final positionX(2) is determined, which is used to calculate the coordinatesr andθ, see 4.30(c). Even though the number of map-pings is determined starting from the leading segment (second point of the crack), the successive mappings start from the last segment which needs to be mapped.
This differs from the approach proposed in [Belytschko and Black, 1999] where the mappings start from leading segment.
A different approach, only of classC0, is to calculate the angleθas θ= arcsind
r, (4.91)
a) Domain under study 1.0
−1.0
−1.5 X1 0.5
X2 Xkink
Xtip
b) MappingC0,αref = 90◦ 1.5
−1.5 1.0
−1.5
−1.0 0.5
X1 X2
F1
c) MappingC1,αref= 90◦ 1.5
−1.5 1.0
−1.5
−1.0
0.5 X1 X2
F1
d) MappingC1,αref = 112.5◦ 1.5
−1.5 1.0
−1.5
−1.0
0.5 X1 X2
F1
Figure 4.27: Comparison of different mappings.
−1.5 −0.5 0.5
−1.5
−1.0 F1
X1
−0.5
C0, αref= 90◦ C1, αref= 90◦ C1, αref= 112.5◦
Figure 4.28: Mappings atX2=−1.0.
where r is calculated directly and dis the signed distance fromX to the crack,
C0, αref= 90◦ C1, αref= 90◦ C1, αref= 112.5◦
−1.5 −0.5 0.5
0.2 0.5 0.9
∂F1
∂X1
X1 Figure 4.29: Partial derivative ∂F1
∂X1 for different mappings atX2=−1.0.
considering that the crack extends from its tips, see for example [Unger et al., 2007].
a)
b)
c) e2
e1
etip1 etip2 X
X(1)
Xkink(1) αref(1) β(1) β(1)rot α(1)rot
α(1)
α(2)
αrot(2)
β(2)
βrot(2)
etip1 etip2 αref(2)
X X(1)
Xtip Xkink(2)
Xtip Xkink(2)
αrot(1)
α(2)rot
X X(1) X(2)
β(1)rot
βrot(2) etip2
etip1
r θ
Figure 4.30: Successive mappings.