Stress Intensity Factor Minimization

The objective of the research was to explore how topology optimization can be used to optimized for damage tolerance objectives such as fatigue crack growth rate. It was hypothesized that the difficulties would lay in the formulation an objective function and the adjoint equation. There formulation should be based upon linear fracture mechanics with the use of FEA.

Continuum formulation

The problem formulation, required for optimization problems, should contain the optimization objective, its link to the design variables and the constraints. Because the goal is design a geometry with the lowest crack growth rate and the Paris-Erdogan law 1 minimizing stress intensity factor \(K_I\) was chosen as the objective. Due to this formulation the design geometry, topology, is the optimization variable.

Fig. 8 Design domain \(\Omega\) with a crack, arbitrary boundary conditions and a density \(X\) which is dependent on the position vector \(\boldsymbol{x}\).

Assuming a general problem, shown in Fig. 8, which minimizes the stress intensity by changing the material distribution, \(X(\boldsymbol{x})\) within the design domain \(\Omega\), the following mathematical formulation is proposed,

\[\begin{split}\min_{X(\boldsymbol{x})} \;\;& K_I(X(\boldsymbol{x}))\\ &\begin{array}{llll} \text{s.t. :} & a(\boldsymbol{u}(X(\boldsymbol{x})),\hat{\boldsymbol{u}}) = l(\hat{\boldsymbol{u}}) \\ & \displaystyle\int_{\Omega} X(\boldsymbol{x}) \text{ d}\Omega \; = \; \text{ Vol}(\Omega^m) \; \leq \; V \\ & X_{\min} \leq X(\boldsymbol{x}) \leq X_{\max} \end{array}\end{split}\]

it enforces equilibrium with a virtual work method while the problem is subjected to a resource constraint. This constraint limits the volume within the design domain that can be filled with a material beside setting a minimum and maximum density value.

For any optimization a link between the objective and the design variables must be made. The method proposed here can be used for two cases, variable thickness plate and discrete material distribution. The honeycomb infill problem is a type of discrete material distribution and will not be discussed separately. In the first case the optimization variables \(X\) are interpreted as the local plate thickness. As the thickness influences the local stiffness properties it affects the stress intensity values at the crack tip. For this variable thickness sheet a linear relation,

\[\boldsymbol{E}_{ijkl}(\boldsymbol{x}) = \boldsymbol{E}_{ijkl, \min} + X(\boldsymbol{x})\left(\boldsymbol{\overline{E}}_{ijkl} - \boldsymbol{E}_{ijkl, \min}\right)\]

between local stiffness and thickness is used. This equation was proposed by M.P. Rossow and J.E. Taylor 2 and discussed by O. Sigmund 3, and causes the stiffness to become twice as high when the thickness is doubled. Here \(\boldsymbol{\overline{E}}_{ijkl}\) is a constant stiffness tensor related to the material it unity thickness while \(\boldsymbol{E}_{ijkl, \min}\) a tensor is with very small stiffness. Which enforces the total stiffness to be larger than zero. One cannot allow the stiffness to become zero as it would cause the FEA to fail. This relation might be inaccurate due to out of plane effects at thickness changes and it will be necessary to measure under what circumstances this equation is invalid.

When the goal is to obtain a discrete design the density values can be either \(0\) (no material) or \(1\) (material). This however causes the objective equation to become discrete as well as the method used a gradient approach and requires a continuous function of density. To ensure a discrete final design while maintaining a continuous objective function a penalization method was implemented. The method used was based upon the penalized proportional stiffness method (SIMP),

\[\boldsymbol{E}_{ijkl}(\boldsymbol{x}) = \boldsymbol{E}_{ijkl, \min} + X(\boldsymbol{x})^p\left(\boldsymbol{\overline{E}}_{ijkl} - \boldsymbol{E}_{ijkl, \min}\right)\]

it causes designs to converge to a \(0\)-\(1\) solution when the penalty factor \(p\) is chosen sufficiently high. Values of \(p\geq 3\) are required for designs to become discrete.

Discretisation

The previous section linked the design variables to the stiffness distribution no official formulation of the stress intensity factors in terms of design variables was made. This formulation is indirectly made through the equilibrium constraint as stiffness distribution influences the stress/displacement field of the loaded part, these stress/displacement distribution can be related to the stress intensity factor. The original equilibrium equation is in a continuum formulation but to simplify the problem a discretized version will be solved using FEA.

To ensure a direct and efficient calculation of the stress intensity factor while using a finite element analysis an enrichment method was used for elements close to the crack tip. The method used was developed by S.E. Benzley 4 and improved by L.N. Gifford 5. It uses a linear summation of a continuous displacement field and a near crack tip displacement field capturing both the discrete behavior at the crack tip and the continuous one around it. The discrete solution was derived with the Westergaard function method 6. This type of tip element enrichment allows accurate predictions of stress intensity directly from the FEA without any post processing as it can be found in the displacement vector.

Crack tip element

The method uses special elements around the crack tip of which the stiffness matrix needs to be derived. As these enriched elements based upon an addition of the continuous and singularity displacement field these are discussed separately at first.

Fig. 9 Nodal definition of the crack tip element.

The enrichment method shown here was based upon the crack tip element developed my L.N. Gifford 5. Who based the enriched elements on a bicubic serendipity elements, see Fig. 9. The algorithm presented here keeps the local coordinate system \((\xi,\, \eta)\) as only a regular mesh with square elements will be used. For a more general element that can contain cracks under an angle and that transforms elements from \((\xi,\, \eta)\) to \((x,\, y)\) see the original paper 5.

The displacement field within the bicubic serendipity 12-node element can be described by:

\[\boldsymbol{u} = \sum_{i=0}^{11} N^i(\xi,\, \eta)\boldsymbol{u}^i\]

where the shape functions \(N^i\) are,

\[\begin{split}&N^0 = \frac{1}{32}\left(1 - \eta\right) \left(1 - \xi\right) \left(9 \eta^{2} + 9 \xi^{2} - 10\right)\\ &N^1 = \frac{9}{32}\left(1 - \eta\right) \left(1 - 3 \xi\right) \left(1 - \xi^{2}\right) \\ &N^2 = \frac{9}{32}\left(1 - \eta\right) \left(1 + 3 \xi\right) \left(1 - \xi^{2}\right)\\ &N^3 = \frac{1}{32}\left(1 - \eta\right) \left(1 + \xi\right) \left(9 \eta^{2} + 9 \xi^{2} - 10\right) \\ &N^4 = \frac{9}{32}\left(1 - 3 \eta\right) \left(1 + \xi\right) \left(1 - \eta^{2}\right) \\ &N^5 = \frac{9}{32}\left(1 + 3 \eta\right) \left(1 + \xi\right) \left(1 - \eta^{2}\right) \\ &N^6 = \frac{1}{32}\left(1 + \eta\right) \left(1 + \xi\right) \left(9 \eta^{2} + 9 \xi^{2} - 10\right)\\ &N^7 = \frac{9}{32}\left(1 + \eta\right) \left(1 + 3 \xi\right) \left(1 - \xi^{2}\right) \\ &N^8 = \frac{9}{32}\left(1 + \eta\right) \left(1 - 3 \xi\right) \left(1 - \xi^{2}\right) \\ &N^9 = \frac{1}{32}\left(1 + \eta\right) \left(1 - \xi\right) \left(9 \eta^{2} + 9 \xi^{2} - 10\right)\\ &N^{10} = \frac{9}{32}\left(1 + 3 \eta\right) \left(1 - \xi\right) \left(1 - \eta^{2}\right)\\ &N^{11} = \frac{9}{32}\left(1 - 3 \eta\right) \left(1 - \xi\right) \left(1 - \eta^{2}\right)\end{split}\]

Added to this will be the crack tip singularity displacement field which derivation starts from the definition of stress intensity factors in a simplified 2D space,

\[\begin{split}K_I = \lim\limits_{r \rightarrow 0} \sqrt{2\pi r} \sigma_{xx}\\ K_{II} = \lim\limits_{r \rightarrow 0} \sqrt{2\pi r} \sigma_{xy}\end{split}\]

and the crack tip stresses derived with the Westergaard method 6,

\[\begin{split}\sigma_{xx} = & \frac{K_I}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\left( 1 - \sin\frac{\theta}{2}\sin\frac{3\theta}{2}\right) \\ &- \frac{K_{II}}{\sqrt{2\pi r}}\sin\frac{\theta}{2}\left(2 + \cos\frac{\theta}{2}\cos\frac{3\theta}{2}\right)\\ \sigma_{yy} =& \frac{K_I}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\left( 1 + \sin\frac{\theta}{2}\sin\frac{3\theta}{2}\right) \\ &+ \frac{K_{II}}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\sin\frac{\theta}{2}\cos\frac{3\theta}{2}\\ \tau_{xy} = & \frac{K_I}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\sin\frac{\theta}{2}\cos\frac{3\theta}{2} \\ &+ \frac{K_{II}}{\sqrt{2\pi r}} \cos\frac{\theta}{2}\left(1 - \sin\frac{\theta}{2}\sin\frac{3\theta}{2}\right)\end{split}\]

which are accurate approximations of the stresses close to the crack tip, i.e. \(r\) is small. Fig. 10 shows the axis system definition for the calculation around the crack tip.

Fig. 10 Definition of the axis systems around the crack tip.

A formulation of the displacement field can be found by integration leading to,

\[\begin{split}u_x =& K_I f_x(r,\, \theta) + K_{II} g_x(r,\, \theta) \\ =& \frac{K_I}{4G}\sqrt{\frac{r}{2\pi}} \left(-1 + \gamma -2\sin^2\frac{\theta}{2}\right)\cos\frac{\theta}{2}\\ &+ \frac{K_{II}}{4G}\sqrt{\frac{r}{2\pi}}\left(1 + \gamma + 2\cos^2\frac{\theta}{2}\right)\sin\frac{\theta}{2}\\ u_y =& K_I f_y(r,\, \theta) + K_{II} g_y(r,\, \theta) \\ =& \frac{K_I}{4G}\sqrt{\frac{r}{2\pi}}\left(1 + \gamma +2\cos^2\frac{\theta}{2} \right)\sin\frac{\theta}{2} \\ &+ \frac{K_{II}}{4G}\sqrt{\frac{r}{2\pi}}\left(1 - \gamma +2\sin^2\frac{\theta}{2}\right)\cos\frac{\theta}{2}\end{split}\]

where \(\gamma = (3-\nu)/(1+\nu)\) for plane stress and \(\gamma = 3-4\nu\) for plane strain 7. When assuming linear fracture mechanics one can describe the displacement field of this element as summation of the continuums and the singularity displacement fields resulting in:

\[\begin{split}u_x = K_I f_x(r,\, \theta) + K_{II}g_x(r,\, \theta) + \sum N^i(\xi,\, \eta)u_x^i \\ u_y = K_I f_y(r,\, \theta) + K_{II}g_y(r,\, \theta) + \sum N^i(\xi,\, \eta)u_y^i\end{split}\]

The singularity equations need to be transformed from the \((r,\, \theta)\) axis into the local \((\xi,\, \eta)\) system. This transformation is dependent of the relative location of the crack tip to the local element axis system.

The enriched displacement functions can cause discontinuities at the border to normal elements, this can be repaired by multiplying the enrichment terms of the displacement function with an equation that is 1 at the crack tip and 0 at the border to non enriched elements 4. It has however been reported that the effects of discontinuities are minor and this solution was therefore not implemented 5.

Following a definition of FE by Zienkiewicz 8 an element stiffness matrix can be calculated with,

\[\boldsymbol{K} = \int_{-1}^{1}\int_{-1}^{1} \boldsymbol{B}^T\boldsymbol{DB}\; \det\boldsymbol{J} \;\; \text{d}\xi \text{d}\eta\]

where \(\boldsymbol{D}\) the material stiffness matrix is, \(\boldsymbol{J}\) the Jacobian of axis system transformation \((\xi,\, \eta)\) into the global \((x,\, y)\) axis system is and \(\boldsymbol{B}\) the matrix is that converts displacement into strain. The integration was performed with a Gauss-Legendre quadrature function with 8x8 integration points as was found sufficient by L.N. Gifford 5.

For a standard bicubic serendipity element this \(\boldsymbol{B}\) matrix is of shape \((3,\, 24)\) however due to the enrichment it becomes \((3,\, 26)\). Which results in a final stiffness matrix of \((26,\, 26)\). Where

\[\begin{split}\boldsymbol{f} = \boldsymbol{K}\boldsymbol{u} = \begin{pmatrix} f_x^0\\ \vdots\\ f^*_x\\ f^*_y \end{pmatrix} = \begin{bmatrix} \boldsymbol{k} & \vdots & \boldsymbol{k}_{12}\\ \dotsm & \vdots & \dotsm \\ \boldsymbol{k}_{21} & \vdots & \boldsymbol{k}_{22} \end{bmatrix} \begin{pmatrix} u_x^0\\ \vdots\\ K_I\\ K_{II} \end{pmatrix}\end{split}\]

Here \(\boldsymbol{k}\) is similar to the stiffness matrix of a normal bicubic element, the enrichment is in the parts \(\boldsymbol{k}_{12}\), \(\boldsymbol{k}_{21}\) and \(\boldsymbol{k}_{22}\). New terms do also appear in the force vector, where \(f^*_x\) and \(f^*_y\) are so-called singular loads. They describe the external forces applied on the crack boundary 4, in general these values are zero.

Meshing strategy

To reduce computational costs these enriched elements are only used at the crack tip and conventional linear elements are used throughout the rest of the mesh. It uses the hanging node method to connect the elements as can be seen in Fig. 11.

Fig. 11 Top section of mesh around a crack tip, \(\oplus\) is the enrichment node with \(K_I\) and \(K_{II}\), while solid circles represent the linear ones and the open circle the higher order ones.

This mesh is not conform which can potentially cause the displacement field to become discontinuous. To avoid this one could use normal bicubic serendipity elements throughout the entire mesh which is computational inefficient. However, using a multi-resolution interpretation of topology optimization its performance might be improved 9.

Currently the linear system of the FEA, \(\boldsymbol{f} = \boldsymbol{Ku}\), and the adjoint equation, \(\boldsymbol{l} = \boldsymbol{K\lambda}\), are solved with a complete Cholesky decomposition. A more efficient methods can be formulated with a Multi Grid Conjugate Gradient method as proposed by O. Amir 10.

Objective formulation

As a spacial discretized method (FEA) was used to calculate the objective the problem formulation needs to become discretized as well. For a mesh of \(N\) elements the optimization objective becomes;

\[\begin{split}&\min_{X_1, X_2, \dots, X_N} \;\; K_I=\boldsymbol{l}^T\boldsymbol{u}\\ &\hspace{0.75cm}\begin{array}{llll} \text{s.t. :} & \boldsymbol{Ku} = \boldsymbol{f} \\ & \displaystyle\sum^N_{e=1} v_eX_e \; \leq \; V \\ & X_{\min} \leq X_e \leq X_{\max} \;\; \forall \;\; e \in \{1, 2, \dots, N\}\\ \text{where :} & \boldsymbol{K} = \displaystyle\sum_{e=1}^{N}\boldsymbol{K}_e(X_e, \overline{E}) \end{array}\end{split}\]

which minimizes the stress intensity factor while ensuring equilibrium and setting constraints to the density distribution. Here \(\boldsymbol{u}\) is the enriched displacement vector, \(f\) the force vector and \(v_e\) is the (relative) element volume. \(\boldsymbol{l}\) is zero vector except for the degree of freedom linked to the stress intensity factor, and the multiplication of \(\boldsymbol{l}^T\boldsymbol{u}\) will return the stress intensity factor. This is similar to the compliant mechanism optimization mentioned by O. Sigmund 11 where the displacement of a specific degree of freedom is maximized.

Sensitivity analysis

The local convex approximation requires the calculation of the sensitivity of \(K_I\) to a density change in any element. This can be measured by \(\partial K_I / \partial X_e\), which can be calculated with the following steps and starts with adding a zero term after the known function \(K_I = \boldsymbol{l}^T\boldsymbol{u}\), where \(\boldsymbol{\lambda}\) is an arbitrary vector:

\[K_I = \boldsymbol{l}^T\boldsymbol{u} - \boldsymbol{\lambda}^T\left(\boldsymbol{Ku} - \boldsymbol{f}\right)\]

\[\frac{\partial K_I}{\partial X_e} = \left(\boldsymbol{l}^T-\boldsymbol{\lambda}^T\boldsymbol{K} \right)\frac{\partial \boldsymbol{u}}{\partial X_e} - \boldsymbol{\lambda}^T\frac{\partial \boldsymbol{K}}{\partial X_e}\boldsymbol{u}\]

Now choosing a convenient vector for \(\boldsymbol{\lambda}\) which causes \(\boldsymbol{l}^T-\boldsymbol{\lambda}^T\boldsymbol{K}\) to be zero leads to the following expression for the sensitivity,

\[\begin{split}\frac{\partial K_I}{\partial X_e} =& - \boldsymbol{\lambda}^T\frac{\partial \boldsymbol{K}}{\partial X_e}\boldsymbol{u}&\\ & \text{where:} \hspace{1cm} \boldsymbol{l} = \boldsymbol{K\lambda}\end{split}\]

This means that \(\boldsymbol{\lambda}\) can be calculated with the FEA, where \(\boldsymbol{l}\) is seen as a sort force vector, by solving \(\boldsymbol{l} = \boldsymbol{Ku}\). The sensitivity of \(\boldsymbol{K}\) to the element density can be calculated, resulting in the following gradient:

\[\frac{\partial K_I}{\partial X_e} = - pX_e^{p-1}\boldsymbol{\lambda}^T\boldsymbol{K}_e\boldsymbol{u}\]

Computational implementation

The iterative implementation of topology optimization as proposed by M. Beckers, 8 or M.P. Bendsøe and O. Sigmund 2 are similar. It exists out of three parts, initialization, optimization and post processing. The flowchart of the local compliance algorithm can be found in Fig. 12.

Fig. 12 Flowchart for fatigue crack growth rate minimization 7.

In the initialization phase the problem is set up. It defines the design domain, the loading conditions, the initial design and generates the finite element mesh that will be used in the optimization phase.

The optimization phase is the iterative method that solves the topology problem. It will analyze the current design with a FEA. After which it will calculate the sensitivity of the stress intensity factor to the density of each element, this is the local gradient of which the calculation is discussed in Sensitivity analysis and MMA. The Method of Moving Asymptotes (MMA), developed by K. Svanberg 9, is used to formulate a simplified convex approximation of the problem which is optimized to formulate the updated design. These steps are performed in a loop until the design is converged, i.e. when the change in design between two iterations becomes negligible.

Post processing is required to remove the last elements with intermediate values and generate a shape out of the design, for example a CAD or STL file. This algorithm will not contain any of the post processing steps. The code used in this communication simply plots the final shape and load case.

Examples and Results

References

1

16. 3. Paris and F. Erdogan, “A Critical Analysis of Crack Propagation Laws,” J. Basic Eng., vol. 85, no. 4, p. 528, 1963.

2(1,2)

13. 16. Rossow and J. E. Taylor, “A Finite Element Method for the Optimal Design of Variable Thickness Sheets,” AIAA J., vol. 11, no. 11, pp. 1566–1569, Nov. 1973.

3

15. Sigmund, N. Aage, and E. Andreassen, “On the (non-)optimality of Michell structures,” Struct. Multidiscip. Optim., vol. 54, no. 2, pp. 361–373, 2016.

4(1,2,3)

19. 5. Benzley, “Representation of singularities with isoparametric finite elements,” Int. J. Numer. Methods Eng., vol. 8, no. 3, pp. 537–545, 1974.

5(1,2,3,4,5)

12. Nash Gifford and P. D. Hilton, “Stress intensity factors by enriched finite elements,” Eng. Fract. Mech., vol. 10, no. 3, pp. 485–496, Jan. 1978.

6(1,2)

8. 13. Westergaard, “Bearing pressures and cracks,” J. Appl. Mech., vol. 6, pp. A49-53, 1939.

7(1,2)

1. 6. Bower, “Modeling Material Failure,” in Applied Mechanics of Solids, 1st ed., Baton Rouge (LA): CRC Press, 2009, pp. 569.

8(1,2)

15. 3. Zienkiewicz, The Finite Element Method In Engineering Science. New York (NY): McGraw-Hill, 1971.

9(1,2)

10. 16. Groen, M. Langelaar, O. Sigmund, and M. Ruess, “Higher-order multi-resolution topology optimization using the finite cell method,” Int. J. Numer. Methods Eng., vol. 110, no. 10, pp. 903–920, Jun. 2017.

10

15. Amir, N. Aage, and B. S. Lazarov, “On multigrid-CG for efficient topology optimization,” Struct. Multidiscip. Optim., vol. 49, no. 5, pp. 815–829, May 2014.

11

15. Sigmund, “On the design of compliant mechanisms using topology optimization,” Mech. Struct. Mach., vol. 25, no. 4, pp. 493–524, 1997.