User:Tpm92/sandbox

Computational Methods In Fracture Mechanics
A handful of methods exist for calculating G with finite elements. Although a direct calculation of the J-integral is possible (using the strains and stresses outputted by FEA), approximate approaches for some type of crack growth exist and provide reasonable accuracy with straightforward calculations. This section will elaborate on some relatively simple methods for fracture analysis utilizing numerical simulations.

Nodal Release Method
[This was done by Christophe]

Includes Figure #1

Modified Crack Closure Integral
Similar to the Nodal Release Method, the Modified Crack Closure Integral (MCCI) is a method for calculating the strain energy release rate utilizing FEA nodal displacements $$(u_i^{j})$$ and forces $$(F_i^{j})$$. Where $$i$$ represents the direction corresponding to the Cartesian basis vectors with origin at the crack tip, and $$j$$ represents the nodal index. MCCI is more computationally efficient than the nodal release method because it only requires one analysis for each increment of crack growth.

A necessary condition for the MCCI method is uniform element length $$(\Delta a )$$ along the crack face in the $$x_1 - $$direction. Additionally, this method requires sufficient discretization such that over the length of one element stress fields are self-similar. This implies that $$K(a+\Delta a) \approx K(a)$$as the crack propagates. Below are examples of the MCCI method with two types of common finite elements.

4-Node Elements
The 4-node square linear elements seen in Figure 2 have a distance between nodes $$j$$ and $$j+1$$ equal to $$\Delta a. $$ Consider a crack with its tip located at node $$j.$$ Similar to the nodal release method, if the crack were to propagate one element length along the line of symmetry (parallel to the $$x_1$$-axis) the crack opening displacement would be the displacement at the previous crack tip, i.e. $$\boldsymbol {u^j}$$and the force at the new crack tip $$(j+1) $$ would be $$\boldsymbol F^{j+1}.$$ Since the crack growth is assumed to be self-similar the displacement at node $$j$$ after the crack propagates is equal to the displacement at node $$j-1$$ before the crack propagates. This same concept can be applied to the forces at node $$j+1$$ and $$j.$$ Utilizing the compliance method as shown in the nodal release section we recover the following equations for strain energy release rate: $$G_1^{\text{MCCI}}=\frac{1}{2\Delta a}F_2^{j}{\Delta u_2^{j-1}} $$

$$G_2^{\text{MCCI}}=\frac{1}{2\Delta a}F_1^{j}{\Delta u_1^{j-1}} $$

$$G_3^{\text{MCCI}}=\frac{1}{2\Delta a}F_3^{j}{\Delta u_3^{j-1}} $$ Where $$\Delta u_i^{j-1} = u_i^{(+)j-1}-u_i^{(-)j-1} $$(displacement above and below the crack face respectively). Because we have a line of symmetry parallel to the crack, we can assume $$u_i^{(+)j-1} = -u_i^{(-)j-1}. $$

Thus, $$\Delta u_i^{j-1} = 2u_i^{(+)j-1}. $$

8-Node Elements
The 8-node rectangular elements seen in Figure 3 have quadratic basis functions. The process for calculating G is the same as the 4-node elements with the exception that $$\Delta a$$ (the crack growth over one element) is now the distance from node $$j$$ to $$j+2.$$ Once again, making the assumption of self-similar straight crack growth the strain energy release rate can be calculated utilizing the compliance method:   $$G_1^{\text{MCCI}}=\frac{1}{2\Delta a}(F_2^{j}{\Delta u_2^{j-2}} + F_2^{j+1}{\Delta u_2^{j-1}})$$

$$G_2^{\text{MCCI}}=\frac{1}{2\Delta a}(F_1^{j}{\Delta u_1^{j-2}} + F_1^{j+1}{\Delta u_1^{j-1}})$$

$$G_3^{\text{MCCI}}=\frac{1}{2\Delta a}(F_3^{j}{\Delta u_3^{j-2}} + F_3^{j+1}{\Delta u_3^{j-1}})$$ Like with the nodal release method the accuracy of MCCI is highly dependent on the level of discretization along the crack tip, i.e. $$G=\lim_{\Delta a \to 0}G^{\text{MCCI}}.$$ Accuracy also depends on element choice. A mesh of 8-node quadratic elements can produce more accurate results than a mesh of 4-node linear elements with the same number of degrees of freedom in the mesh.

Domain Integral Approach for J
"[This was done by Christophe]"Includes Figure #4

2-D Crack Tip Singular Elements
The above-mentioned methods for calculating strain energy release rate asymptotically approach the actual solution with increased discretization, but fail to fully capture the crack tip singularity. More accurate simulations can be preformed by utilizing quarter-point elements around the crack tip. These elements have a built in singularity which more accurately produces stress fields around the crack tip. The advantage of the quarter-point method is that it allows for coarser finite element meshes and greatly reduces computational cost. Furthermore, these elements are derived from small modifications to common finite elements without requiring special computational programs for analysis. For the purposes of this section elastic materials will be examined, although this method can be extended to elastic-plastic fracture mechanics. Assuming perfect elasticity the stress fields will experience a $$\frac{1}{\sqrt r} $$ crack tip singularity.

8-Node Isoparametric Element
The 8-node quadratic element are described by Figure 5 in both parent space with local coordinates $$\xi$$ and $$\eta,$$ and by the mapped element in physical/global space by $$x$$ and $$y.$$ The parent element is mapped from the local space to the physical space by the shape functions $$N_i(\xi,\eta)$$ and the degree of freedom coordinates $$(x_i,y_i). $$ The crack tip is located at $$\xi = -1,\eta = -1 $$ or $$x = 0,y = 0. $$ $$x(\xi,\eta) = \sum_{i=1}^{8}N_i(\xi,\eta)x_i $$  $$y(\xi,\eta) = \sum_{i=1}^{8}N_i(\xi,\eta)y_i $$

In a similar way, displacements (defined as $$u\equiv u_1,v\equiv u_2 $$) can also be mapped. $$u(\xi,\eta) = \sum_{i=1}^8 N_i(\xi,\eta)u_i $$ $$v(\xi,\eta) = \sum_{i=1}^8 N_i(\xi,\eta)v_i $$ A property of shape functions in the finite element method is compact support, specifically the Kronecker delta property (i.e. $$N_i = 1 $$ at node $$i $$ and zero at all other nodes). This results in the following shape functions for the 8-node quadratic elements : $$N_1 = \frac{-(\xi-1)(\eta-1)(1+\eta+\xi)}{4} $$

$$N_2 = \frac{(\xi+1)(\eta-1)(1+\eta-\xi)}{4} $$

$$N_3 = \frac{(\xi+1)(\eta+1)(-1+\eta+\xi)}{4} $$

$$N_4 = \frac{-(\xi-1)(\eta+1)(-1+\eta-\xi)}{4} $$

$$N_5 = \frac{(1-\xi^2)(1-\eta)}{2} $$

$$N_6 = \frac{(1+\xi)(1-\eta^2)}{2} $$

$$N_7 = \frac{(1-\xi^2)(1+\eta)}{2} $$

$$N_8 = \frac{(1-\xi)(1-\eta^2)}{2} $$ When considering a line in front of the crack that is co-linear with the $$x$$- axis (i.e. $$N_i(\xi,\eta=-1)$$) all basis functions are zero except for $$N_{1,2,5}.$$ $$N_1(\xi,-1) = -\frac{\xi(1-\xi)}{2} $$

$$N_2(\xi,-1) = \frac{\xi(1+\xi)}{2} $$

$$N_5(\xi,-1) = (1-\xi^2) $$ Calculating the normal strain involves using the chain rule to take the derivative of displacement with respect to $$x. $$ $$\gamma_{xx} = \frac{\partial u}{\partial x} = \sum_{i = 1,2,5}\frac{\partial N_i}{\partial \xi}\frac{\partial \xi}{\partial x}u_i $$

If the nodes are spaced evenly on the rectangular element then the strain will not contain the singularity. By moving nodes 5 and 8 position to a quarter of the length $$(\tfrac{L}{4}) $$ of the element closer to the crack tip as seen in figure 5, the mapping from $$\xi \rightarrow x $$ becomes: $$x(\xi) = \frac{\xi(1+\xi)}{2}L + (1-\xi^2)\frac{L}{4} $$ Solving for $$\xi $$ and taking the derivative results in: $$\xi(x) = -1+2\sqrt\frac{x}{L} $$ $$\frac{\partial\xi}{\partial x} = \frac{1}{\sqrt{xL}} $$ Plugging this result into the equation for strain the final result is obtained: $$\gamma_{xx} = \frac{4}{L}(\frac{u_2}{2} - u_5) + \frac{1}{\sqrt{xL}}(2u_5-\frac{u_2}{5}) $$ By moving the mid-nodes to a quarter position results in the correct $$\frac{1}{\sqrt r} $$ crack tip singularity.

Other Element Types
The rectangular element method does not allow for singular elements to be easily meshed around the crack tip. This impedes the ability to capture the angular dependence of the stress fields which is critical in determining the crack path. Also, except along the element edges the $$\frac{1}{\sqrt{r}}$$ singularity exists in a very small region near the crack tip. Figure 6 shows another quarter-point method for modeling this singularity. The 8-node rectangular element can be mapped into a triangle. This is done by collapsing the nodes on the line $$\xi = -1$$ to the mid-node location and shifting the mid-nodes on $$\eta = \pm1$$ to the quarter-point location. The collapsed rectangle can more easily surround the crack tip but requires that the element edges be straight or the accuracy of calculating the stress intensity factor will be reduced.

A better candidate for the quarter-point method is the natural triangle as seen in Figure 7. The element's geometry allows for the crack tip to be easily surrounded and meshing is simplified. Following the same procedure described above, the displacement and strain field for the triangular elements are: "$u = u_3+\sqrt{\frac{x}{L}}[4u_6-3u_3-u_1] + \frac{x}{L}[2u_1+2u_3-4u_6]$|undefined" $$\gamma_{xx} = \frac{\partial u}{\partial x} = \frac{1}{\sqrt{xL}}[-\frac{u_1}{2}-\frac{3u_3}{2}+2u_6]+ \frac{1}{L}[2u_1+2u_3-4u_6]$$ This method reproduces the first two terms of the Williams solutions with a constant and singular term.

An advantage of quarter-point method is that it can be easily generalized to 3-dimensional models. This can greatly reduce computation when compared to other 3-dimensional methods but can lead to error if that crack tip propagates with a large degree of curvature.