User:Cbonneville/sandbox

Computational Methods In Fracture Mechanics

Many 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 finite element analysis), approximate approaches for some types of crack growth exist and provide reasonable accuracy with straightforward calculations.

Nodal Release Method

Figure 1: Crack debounding between two consecutive time steps using nodal release (NR)

If the crack is growing straight, the energy release rate can be decomposed as a sum of 3 terms $G_{i}$associated with the energy in each 3 modes. As a result, the Nodal Release method (NR) can be used to determine $G_{i}$from FEA results. The energy release rate is calculated at the nodes of the finite element mesh for the crack at an initial length and extended by a small distance $\Delta a$. First, we calculate the displacement variation at the node of interest $\Delta {\vec {u}}={\vec {u}}^{(t+1)}-{\vec {u}}^{(t)}$(before and after the crack tip node is released). Secondly, we keep track of the nodal force ${\vec {F}}$ outputted by FEA. Finally, we can find each components of $G$ using the following formulas:

$$G_{1}^{\text{NR}}={\frac {1}{\Delta a}}F_{2}{\frac {\Delta u_{2}}{2}}$$$$G_{2}^{\text{NR}}={\frac {1}{\Delta a}}F_{1}{\frac {\Delta u_{1}}{2}}$$$$G_{3}^{\text{NR}}={\frac {1}{\Delta a}}F_{3}{\frac {\Delta u_{3}}{2}}$$Where $\Delta a$ is the width of the element bounding the crack tip. The accuracy of the method highly depends on the mesh refinement, both because the displacement and forces depend on it, and because $G=\lim _{\Delta a\to 0}G^{\text{NR}}$. Note that the equations above are derived using the crack closure integral.

If the energy release rate exceed a critical value, the crack will grow. In this case, a new FEA simulation is perfomed (for the next time step) where the node at the crack tip is released. For a bounded substrate, we may simply stop enforcing fixed Dirichlet boundary conditions at the crack tip node of the previous time step (i.e. displacements are no longer restrained) . For a symmetric crack, we would need to update the geometry of the domain with a longer crack opening (and therefore generate a new mesh).

MCCI Method

Terrence is doing this part

Domain Integral Approach for J

Contour domain integral for the J-Integral

The J-integral may be calculated directly using the finite element mesh and shape functions. We consider a domain countour as shown in figure 2 and choose an arbitrary smooth function ${\tilde {q}}(x_{1},x_{2})=\sum _{i}N_{i}(x_{1},x_{2}){\tilde {q}}_{i}$ such that ${\tilde {q}}=1$ on $\Gamma$ and ${\tilde {q}}=0$ on ${\mathcal {C}}_{1}$.

For linear elastic cracks growing straight ahead, $G=J$. The energy release rate can then be calculated over the area bounded by the contour using an updated formulation:

$J=\int _{\mathcal {A}}(\sigma _{ij}u_{i,1}{\tilde {q}}_{,j}-W{\tilde {q}}_{,1})d{\mathcal {A}}$

The formula above may be applied to any annular area surrounding the crack tip (in particular, a set of neighboring elements can be used). This method is very accurate, even with a coarse mesh arround the crack tip (one may choose an integration domain located far away, with stresses and displacement less sensitive to mesh refinment)

Derivation of the J-integral for domain integral approach

The J-intregral may be expressed over the full contour as follow: $G=J=\int _{\Gamma }(Wn_{1}-t_{i}{\frac {\partial u_{i}}{\partial x_{1}}})d\Gamma =\int _{{\mathcal {C}}_{1}+{\mathcal {C}}_{+}+{\mathcal {C}}_{-}-{\mathcal {C}}}(Wn_{1}-\sigma _{ij}n_{j}{\frac {\partial u_{i}}{\partial x_{1}}}){\tilde {q}}d\Gamma$ With ${\mathcal {C}}={\mathcal {C}}_{1}+{\mathcal {C}}_{+}+{\mathcal {C}}_{-}-\Gamma$. ${\vec {n}}=-{\vec {m}}$, ${\tilde {q}}=0$ on ${\mathcal {C}}_{1}$and the work and stresses cancel out on ${\mathcal {C}}_{+}$and ${\mathcal {C}}_{-}$, hence by application of the divergence theorem this leads to: $J=\int _{\mathcal {C}}(-Wm_{1}+\sigma _{ij}m_{i}{\frac {\partial u_{i}}{\partial x_{1}}})d{\mathcal {C}}=\int _{\mathcal {A}}({\frac {\partial \sigma _{ij}u_{i,1}{\tilde {q}}}{\partial x_{j}}}-{\frac {\partial W{\tilde {q}}}{\partial x_{1}}})d{\mathcal {A}}$ Finally, by noting that ${\frac {\partial W}{\partial x_{1}}}=\sigma _{ij}{\frac {\partial ^{2}u_{i}}{\partial x_{j}\partial x_{1}}}$and using the equilibrium equation: $J=\int _{\mathcal {A}}(\sigma _{ij}u_{i,1}{\tilde {q}}_{,j}-W{\tilde {q}}_{,1})d{\mathcal {A}}$

Crack Tip Singular Elements

Terrence is doing this part