Topology optimization for transversely isotropic materials with high-cycle fatigue as a constraint

Abstract

We propose a topology optimization method for design of transversely isotropic elastic continua subject to high-cycle fatigue. The method is applicable to design of additive manufactured components, where transverse isotropy is often manifested in the form of a lower Young’s modulus but a higher fatigue strength in the build direction. The fatigue constraint is based on a continuous-time model in the form of ordinary differential equations governing the time evolution of fatigue damage at each point in the design domain. Such evolution occurs when the stress state lies outside a so-called endurance surface that moves in stress space depending on the current stress and a back-stress tensor. Pointwise bounds on the fatigue damage are approximated using a smooth aggregation function, and the fatigue sensitivities are determined by the adjoint method. Several problems where the objective is to minimize mass are solved numerically. The problems involve non-periodic proportional and non-proportional load histories. Two alloy steels, AISI-SAE 4340 and 34CrMo6, are treated and the respective as well as the combined impact of transversely isotropic elastic and fatigue properties on the design are compared.

1 Introduction

Additive manufacturing (AM) is a versatile manufacturing process in which a structural component is fabricated layer-by-layer from digital information. This form of manufacturing is gaining popularity in aerospace and automotive industries. Through AM one can fabricate highly complex structures and together with topology optimization (TO), a powerful design tool, it shares the property of providing a very large freedom in geometrical shape.

The properties of AM-fabricated components are anisotropic due to the layer-wise generation. For example, it has been observed that Ti6Al4V tensile test specimens built using electron beam melting have superior strength and elastic moduli in flat-built specimens compared with top-built specimens (Ladani et al. 2014), and according to Kumara et al. (2018) the AM material alloy 718 exhibits transversely isotropic elastic properties with the lowest Young’s modulus value in the build direction. An important design criterion that needs to be considered in AM-fabricated parts is the fatigue behaviour. Several important factors such as surface roughness and build orientation affect the design against fatigue in AM (Molaei and Fatemi 2018). For instance, AM components used with as-built surface condition (without any post surface treatments) has significantly worse fatigue-life compared with a post surface finish AM components.

In this work we consider high-cycle fatigue (HCF) for materials with transversely isotropic properties in the TO formulation. We incorporate a continuous-time HCF model, in which the stresses are decomposed into longitudinal and transverse directions (Holopainen et al. 2016). This HCF model can handle arbitrary load histories, including non-proportional loads, without use of any cycle-counting algorithm.

Examples of TO formulations with fatigue constraints include Holmberg et al. (2014), where the fatigue constraints were implemented as stress constraints, and Jeong et al. (2018), Collet et al. (2017), and Oest and Lund (2017) which use cycle-counting algorithms. Other examples of structural optimization with fatigue constraints are found in Oest et al. (2017), where the fatigue prediction is done for a simplified damage model assuming a periodic load, and in Gerzen et al. (2017), where sizing optimization is done with fatigue constraints at the welded joints using a cycle-counting algorithm (rainflow-counting). However, these models are restricted to proportional loading histories. Reference Zhang et al. (2019) uses classical techniques, including rainflow-counting, mean stress correction and the Palmgren-Miner rule. It treats fatigue life induced by non-proportional loads and is implemented in TO problems. Our recent contribution (Suresh et al. 2020) implements an evolution-based fatigue model in TO problems. The fatigue model is capable of handling arbitrary load histories, including non-proportional loads. However, all of the above-mentioned optimization formulations with fatigue constraints are developed for materials with isotropic properties. Hence, an extension of a fatigue model in TO that can handle anisotropic materials properties is needed for design of AM parts. This motivated us to extend our previous contribution (Suresh et al. 2020) and implement a HCF model for transversely isotropic materials in a TO problem; in this case minimization of mass subject to a HCF constraint.

The continuous-time fatigue model in Ottosen et al. (2008) uses a so-called endurance surface that moves in the stress space. The model is based on ordinary differential equations (ODEs) that govern the time evolution of fatigue damage at each point in the design domain. Damage development only occurs if the stress state lies outside the endurance surface. The flexibility of the model has led to further developments such as Brighenti et al. (2011), where the fatigue is assessed for complex multiaxial load histories, and Ottosen et al. (2018), where the multiaxial fatigue criterion considers stress gradient effects in critical regions like holes and notches. The model has also been extended to account for anisotropic properties in Holopainen et al. (2016), in particular transverse isotropy. Furthermore, the validity range and computational acceleration of the continuous-time HCF model are investigated in Lindström et al. (2020).

In the present work, HCF is implemented as a constraint on the maximum fatigue damage in the design domain (approximated by an aggregation function, namely the p-norm (Kennedy and Hicken 2015)). With the HCF model, the geometry and material properties remain unaffected, i.e. the finite element (FE) stiffness matrix is constant throughout the load history for a given design. Considering transversely isotropic material properties, we set the lowest Young’s modulus value in the build direction and evaluate the stress history from an anisotropic elastic analysis. Then these stresses are used in the transversely isotropic HCF model to compute the total accumulated fatigue damage by solving the ODEs of the damage and back-stress.

2 Continuous-time fatigue model

We are concerned with HCF, which means that fatigue damage neither influences the geometry, nor the material properties. Assuming small deformations, the constitutive response is considered as linear elastic. For each design, the stress evolution σ(t), with time t belonging to the time interval [0,T], is computed first, and is then used to estimate the fatigue damage through the continuous-time fatigue model suggested in Ottosen et al. (2008) and Holopainen et al. (2016).

A moving endurance surface {σ | β(σ,α) = 0}, where α(t) is the back-stress tensor and β is the endurance function, is used to formulate this fatigue model. It is assumed that the development of fatigue damage occurs only if the stress state lies outside the endurance surface, i.e. β > 0, and the endurance function value is increasing, i.e. $\dot {\beta } > 0$, where a superposed dot indicates time derivative.

Considering a transversely isotropic material, both the elastic properties and the fatigue sensitivity may exhibit directional dependencies. In AM parts, the lowest value of Young’s modulus is in the build direction (Kumara et al. 2018) and the transversely isotropic material matrix shown in Appendix can be used. The fatigue model in Holopainen et al. (2016) is used to account for the anisotropy of the fatigue properties. The key distinguishing feature of this fatigue model is the way the endurance function β is formulated.

From Ottosen et al. (2008) and Suresh et al. (2020), we take the evolution equations for the back-stress α and fatigue damage D as

$$ \begin{array}{@{}rcl@{}} &\dot{\boldsymbol{\alpha}} = C({\boldsymbol{s}} - {\boldsymbol{\alpha}})H(\beta)H(\dot{\beta})\dot{\beta}, \quad &{\boldsymbol{\alpha}}(0) = {\boldsymbol{0}}, \end{array} $$

(1)

$$ \begin{array}{@{}rcl@{}} &\dot{D} = KH(\beta)H(\dot{\beta})\exp(L\beta)\dot{\beta}, \quad &D(0) = 0, \end{array} $$

(2)

where C > 0,K > 0 and L > 0 are material parameters. The deviatoric stress is ${\boldsymbol {s}} = {\boldsymbol {\sigma }} - \frac {1}{3}\text {tr}({\boldsymbol {\sigma }}){\boldsymbol {I}}$, with I as the second-order identity tensor. H is the Heaviside step function and tr(⋅) is the trace of a tensor. The damage D(t) at a point is a real-valued scalar that increases from D = 0, being undamaged material, to D = 1 as critical failure. It remains to formulate the endurance function β and calculate its rate.

2.1 Transversely isotropic HC fatigue model

The fatigue model in Holopainen et al. (2016) considers transversely isotropic material properties. A unit vector m is introduced in the longitudinal (build) direction and the structural tensor is written as M = m ⊗m, where ⊗ is the tensor product symbol. The stress is additively decomposed into longitudinal σ_∥ and transverse σ_⊥ stress tensors, i.e.

$$ {\boldsymbol{\sigma}} = {\boldsymbol{\sigma}}_{\parallel} + {\boldsymbol{\sigma}}_{\bot}, $$

(3)

where the transverse component σ_⊥ is calculated by the use of a projection tensor P defined as P = I −M. This gives

$$ {\boldsymbol{\sigma}}_{\bot} = {\boldsymbol{P}}{\boldsymbol{\sigma}}{\boldsymbol{P}} = {\boldsymbol{\sigma}} - {\boldsymbol{M}}{\boldsymbol{\sigma}} - {\boldsymbol{\sigma}}{\boldsymbol{M}} + {\boldsymbol{M}}{\boldsymbol{\sigma}}{\boldsymbol{M}},\\ $$

and

$$ {\boldsymbol{\sigma}}_{\parallel} = {\boldsymbol{\sigma}} - {\boldsymbol{\sigma}}_{\bot} = {\boldsymbol{M}}{\boldsymbol{\sigma}} + {\boldsymbol{\sigma}}{\boldsymbol{M}} - {\boldsymbol{M}}{\boldsymbol{\sigma}}{\boldsymbol{M}}. $$

To formulate the endurance function we introduce four tensor invariants,

$$ \begin{aligned} &I_{1} = \text{tr}({\boldsymbol{\sigma}}), \quad &I_{2} = \frac{1}{2}\text{tr}({\boldsymbol{\sigma}}^{2}),\\ &I_{4} = \text{tr}({\boldsymbol{\sigma}}{\boldsymbol{M}}), \quad &I_{5} = \text{tr}({\boldsymbol{\sigma}}^{2}{\boldsymbol{M}}). \end{aligned} $$

Following Holopainen et al. (2016), the endurance function is written as

$$ \begin{aligned} \beta({\boldsymbol{\sigma}},{\boldsymbol{\alpha}}) = \frac{1}{S_{\bot}}\left[\bar{\sigma} + A_{\parallel}\text{tr}({\boldsymbol{\sigma}}_{\parallel}) + A_{\bot}\text{tr}({\boldsymbol{\sigma}}_{\bot}) - \tilde{S}_{0}({\boldsymbol{\sigma}})\right], \end{aligned} $$

(4)

where A_∥ > 0 and A_⊥ > 0 are material parameters related to the corresponding physical directions. It can be shown that tr(σ_∥) = I₄ ¹ and tr(σ_⊥) = I₁ − I₄. The effective stress $\bar {\sigma }$ is obtained by

$$ \begin{aligned} \bar{\sigma} = \sqrt{\frac{3}{2}({\boldsymbol{s}} - {\boldsymbol{\alpha}}):({\boldsymbol{s}} - {\boldsymbol{\alpha}})}, \end{aligned} $$

where the colon operator ‘:’ denotes the Frobenius inner product, i.e. A : D = tr(AD^T), with A and D as second-order tensors. The interpolation of endurance stress is $\tilde {S}_{0}({\boldsymbol {\sigma }}) = (1-\xi ({\boldsymbol {\sigma }}))S_{\bot } + \xi ({\boldsymbol {\sigma }}) S_{\parallel }$, with S_∥ > 0 and S_⊥ > 0 as the endurance stresses in the corresponding directions. The scalar function ξ is the stress ratio, i.e.

$$ \begin{aligned} \xi({\boldsymbol{\sigma}}) = \frac{{\boldsymbol{\sigma}}_{\parallel}:{\boldsymbol{\sigma}}_{\parallel}}{{\boldsymbol{\sigma}}:{\boldsymbol{\sigma}}} = \frac{2I_{5} - {I_{4}^{2}}}{2I_{2}}. \end{aligned} $$

(5)

Taking the time derivative of (4) and with help of (5), we get

$$ \begin{array}{@{}rcl@{}} \begin{aligned} \dot{\beta} = &\frac{1}{S_{\bot}}\left[\frac{3}{2}\frac{({\boldsymbol{s}} - {\boldsymbol{\alpha}})}{\bar{\sigma}} {}:{} \dot{{\boldsymbol{s}}} - \frac{3}{2}\frac{({\boldsymbol{s}} - {\boldsymbol{\alpha}})}{\bar{\sigma}} {}:{} \dot{\boldsymbol{\alpha}} + (A_{\parallel} - A_{\bot}) \boldsymbol{M}:\dot{\boldsymbol{\sigma}}\right. \\+ & A_{\bot}{\boldsymbol{I}}:\dot{{\boldsymbol{\sigma}}} - \frac{(S_{\parallel} - S_{\bot})}{2{I_{2}^{2}}}\left( I_{2}(4{\boldsymbol{\sigma}} - 2I_{4}{\boldsymbol{I}}){\boldsymbol{M}} - \right.\\ &\left.\left. (2I_{5} - {I_{4}^{2}}){\boldsymbol{\sigma}}\right):\dot{{\boldsymbol{\sigma}}}\right]. \end{aligned} \end{array} $$

Substituting (1) for $\dot {\boldsymbol {\alpha }}$, and rearranging, we have

$$ \begin{aligned} \underset{>0}{\underbrace{\left( S_{\bot} + H(\beta)H(\dot{\beta})C\bar{\sigma}\right)}}\dot{\beta} &= \left[\frac{3}{2}\frac{({\boldsymbol{s}} - {\boldsymbol{\alpha}})}{\bar{\sigma}} + (A_{\parallel} - A_{\bot}) \boldsymbol{M} + \right.\\ & A_{\bot}{\boldsymbol{I}} -\frac{(S_{\parallel} - S_{\bot})}{2{I_{2}^{2}}}\left( I_{2}(4{\boldsymbol{\sigma}} - 2I_{4}{\boldsymbol{I}}){\boldsymbol{M}} - \right. \\ & \left.\left. (2I_{5} - {I_{4}^{2}}){\boldsymbol{\sigma}}\right)\right]:\dot{{\boldsymbol{\sigma}}}. \end{aligned} $$

When $H(\dot {\beta }) > 0$ and since $\left (S_{\bot } + H(\beta )H(\dot {\beta })C\bar {\sigma }\right ) > 0$, we can write the above rate function as

$$ \begin{aligned} \dot{\beta}({\boldsymbol{\sigma}},{\boldsymbol{\alpha}},\dot{{\boldsymbol{\sigma}}}) = &\frac{1}{S_{\bot} + H(\beta)C\bar{\sigma}}\left[\frac{3}{2}\frac{({\boldsymbol{s}} - {\boldsymbol{\alpha}})}{\bar{\sigma}} + (A_{\parallel} - A_{\bot}) \boldsymbol{M} + \right.\\ & A_{\bot}{\boldsymbol{I}} -\frac{(S_{\parallel} - S_{\bot})}{2{I_{2}^{2}}}\left( I_{2}(4{\boldsymbol{\sigma}} - 2I_{4}{\boldsymbol{I}}){\boldsymbol{M}} - \right. \\ & \left. \left. (2I_{5} - {I_{4}^{2}}){\boldsymbol{\sigma}}\right)\right]:\dot{{\boldsymbol{\sigma}}}. \end{aligned} $$

(6)

The isotropic HCF model in Ottosen et al. (2008) emerges as a special case when A_⊥ = A_∥ = A > 0 and S_⊥ = S_∥ = S₀ > 0. Then the endurance function (4) and its rate (6) become

$$ \begin{array}{@{}rcl@{}} \beta({\boldsymbol{\sigma}},{\boldsymbol{\alpha}})& =& \frac{1}{S_{0}}\left[\bar{\sigma} + A \text{tr}({\boldsymbol{\sigma}}) - S_{0}\right], \\ \dot{\beta}({\boldsymbol{\sigma}},{\boldsymbol{\alpha}},\dot{{\boldsymbol{\sigma}}}) &=& \frac{1}{S_{0} + H(\beta)C\bar{\sigma}}\left[\frac{3}{2}\frac{({\boldsymbol{s}} - {\boldsymbol{\alpha}})}{\bar{\sigma}} + A{\boldsymbol{I}}\right]:\dot{{\boldsymbol{\sigma}}}. \end{array} $$

With the endurance function and its rate defined in (4) and (6), the fatigue damage is estimated by integrating the ODEs (1) and (2).

2.2 Discretization

The transversely isotropic HC fatigue problem is solved numerically. The time domain [0,T] is divided into a finite number of time steps of equal length Δt, i.e. t_i = iΔt, with i = 0,1,2,...,N. The evolution (1) and (2) are approximated by the forward Euler scheme. The stresses evaluated at different steps are denoted σ_i = σ(iΔt), s|_i = s(σ_i), σ_∥|_i = σ_∥(σ_i) and σ_⊥|_i = σ_⊥(σ_i), where |_i denotes function evaluation at the time step i.

The discrete versions of (1) and (2) are written as

$$ \begin{array}{@{}rcl@{}} &{\boldsymbol{\alpha}}_{i} - {\boldsymbol{\alpha}}_{i-1} = C({\boldsymbol{s}}|_{i-1} - {\boldsymbol{\alpha}}_{i-1})H(\beta|_{i-1}){\Delta}{\beta}|_{i}, \end{array} $$

(7)

$$ \begin{array}{@{}rcl@{}} &{D}_{i} - {D}_{i-1} = K\exp(L\beta |_{i-1})H(\beta |_{i-1}){\Delta}{\beta} |_{i}. \end{array} $$

(8)

The increment of the endurance function in (6) is approximated as

$$ \begin{aligned} {\Delta} \beta |_{i} = &{\Delta} t \dot{\beta}\left( {\boldsymbol{\sigma}}_{i-1}, \frac{{\boldsymbol{\sigma}}_{i}-{\boldsymbol{\sigma}}_{i-1}}{\Delta t}, {\boldsymbol{\alpha}}_{i-1}\right)\\ = &\frac{1}{S_{\bot} + H(\beta |_{i-1})C\bar{\sigma} |_{i-1}}\left[\frac{3}{2}\frac{({\boldsymbol{s}} |_{i-1} - {\boldsymbol{\alpha}}_{i-1})}{\bar{\sigma} |_{i-1}} + A_{\parallel}{\boldsymbol{M}}-\right.\\ &A_{\bot}{\boldsymbol{M}} + A_{\bot}{\boldsymbol{I}} -\frac{(S_{\parallel} - S_{\bot})}{2I_{2} |_{i-1}^{2}}\left( I_{2} |_{i-1}(4{\boldsymbol{\sigma}} |_{i-1} - 2I_{4} |_{i-1}{\boldsymbol{I}}){\boldsymbol{M}} \right.\\ &\left. \left.-(2I_{5} |_{i-1} - I_{4} |_{i-1}^{2}){\boldsymbol{\sigma}} |_{i-1}\right)\right]:{\Delta}{{\boldsymbol{\sigma}}}_{i}, \end{aligned} $$

(9)

with the stress increment as Δσ|_i = σ_i −σ_i− 1, and the effective stress as $\bar {\sigma } |_{i} = \bar {\sigma }({{\boldsymbol {\sigma }}_{i}}, {\boldsymbol {\alpha }}_{i})$. The invariants are computed as I₂|_i = I₂(σ_i), I₄|_i = I₄(σ_i) and I₅|_i = I₅(σ_i). The expression (9) is only evaluated if the Frobenius inner product between the terms in the square parenthesis and the stress increment is positive, otherwise Δβ_i = 0.

3 Optimization problem formulation

Through the FE method we obtain a spatially discretized model in which the design domain Ω (Fig. 1) is divided into n_e finite elements. We use a standard, density-based TO method, where the design variables x_e, e = 1,2,...,n_e, are collected in a vector x. Each design variable x_e is associated to a single FE. These variables should ideally take the value 0 (void) or 1 (material) to get so-called black-and-white designs. A design variable filter (Bruns and Tortorelli 2001) is applied to avoid mesh dependency and checker-board patterns arising for commonly used low-order FEs. This filter gives a new set of variables ${\boldsymbol {\rho }} = {\boldsymbol {\rho }}({\boldsymbol {x}}) = [\rho _{1}({\boldsymbol {x}}), \rho _{2}({\boldsymbol {x}}), ..., \rho _{n_{e}}({\boldsymbol {x}})]^{\text {T}}$ called physical variables since they define the structural stiffness and the mass.

Fig. 1

Loading conditions with two load cases: a Load case 1 contains a static load F, b load case 2 contains a dynamic load history $\tilde {{\boldsymbol {F}}}(t)$

To account for both stiffness and fatigue, we create two load cases, shown in Fig. 1. In the first load case (Fig. 1a) we take a static load F to compute the compliance (inverse of stiffness), while in the second load case (Fig. 1b) a dynamic load history $\tilde {{\boldsymbol {F}}}(t)$ is used for estimating the fatigue damage.

The equilibrium equation for the static load case is

$$ {\boldsymbol{K}}({\boldsymbol{\rho}}({\boldsymbol{x}})){\boldsymbol{u}} = {\boldsymbol{F}}. $$

For the fatigue load case, we assume quasi-static equilibrium, and consider a sequence of static problems,

$$ {\boldsymbol{K}}({\boldsymbol{\rho}}({\boldsymbol{x}}))\tilde{{\boldsymbol{u}}}(t) = \tilde{{\boldsymbol{F}}}(t), \quad \forall\ t\in[0,T]. $$

In the equilibrium equations, K(ρ(x)) is the global stiffness matrix, and u and $\tilde {{\boldsymbol {u}}}(t)$ are the displacement vectors for the corresponding load cases.

At each time step t_i = iΔt, we have

$$ {\boldsymbol{K}}({\boldsymbol{\rho}}({\boldsymbol{x}}))\tilde{{\boldsymbol{u}}}_{i} = \tilde{{\boldsymbol{F}}}_{i}, $$

(10)

where $\tilde {{\boldsymbol {F}}}_{i} = \tilde {{\boldsymbol {F}}}(t_{i})$, and the solution for a given design x and time step i is denoted by $\tilde {{\boldsymbol {u}}}_{i}({\boldsymbol {x}})$.

To promote black-and-white designs, the SIMP approach (Bendsoe and Sigmund 2004; Christensen and Klarbring 2009) is used, i.e.

$$ \begin{aligned} &{\boldsymbol{K}}({\boldsymbol{\rho}}({\boldsymbol{x}})) = \overset{n_{e}}{\underset{e=1}{\sum}}(\rho_{e}({\boldsymbol{x}}))^{q}{\boldsymbol{K}}_{e}, \end{aligned} $$

(11)

where K_e are elemental stiffness matrices and q > 1 is a penalization parameter.

For a given $\tilde {{\boldsymbol {u}}}_{i}({\boldsymbol {x}})$, the stress vector at a stress evaluation point is calculated as

$$ \begin{aligned} &\hat{{\boldsymbol{\sigma}}}_{i}({\boldsymbol{x}}) = \hat{{\boldsymbol{\sigma}}}_{i}({\boldsymbol{x}}, \tilde{{\boldsymbol{u}}}_{i}({\boldsymbol{x}})) = {\boldsymbol{E}}{\boldsymbol{B}}\tilde{{\boldsymbol{u}}}_{i}({\boldsymbol{x}}), \end{aligned} $$

(12)

where E is the constitutive tensor in Voigt form, shown in Appendix, and B is the strain-displacement matrix. The stress used to compute the fatigue response is then given by

$$ \begin{aligned} &{\boldsymbol{\sigma}}_{i}({\boldsymbol{x}}) = \left( \frac{\rho_{e}({\boldsymbol{x}}) - \epsilon}{1 - \epsilon}\right)^{r}\hat{{\boldsymbol{\sigma}}}_{i}({\boldsymbol{x}}), \end{aligned} $$

(13)

where r < q is a parameter introduced to avoid the stress singularity phenomenon (Bruggi 2008; Holmberg et al. 2013), and 0 < 𝜖 ≪ 1 is the lower bound on the design variable. The penalization factor gives exactly zero stress in voids. Hence, no artificial damage is developed in such regions.

We are concerned with mass minimization of the component, where the total structural mass is minimized subject to fatigue constraints and a compliance constraint. The compliance is evaluated from the static load case as F^Tu(x). As for the fatigue load case, for a given load history $\tilde {{\boldsymbol {F}}}$ and assuming a single stress evaluation point in each element, the fatigue damage evaluated at a time step i is D_i,e(x). The fatigue constraint for each element is formulated using the total accumulated damage D_N,e(x), i.e. the fatigue damage evaluated at the final time step T = t_N = NΔt, as

$$ D_{N,e}({\boldsymbol{x}}) \le \bar{D}_{e}, \quad e = 1, 2, ..., n_{e}, $$

(14)

with $\bar {D}_{e}$ as the maximum allowable fatigue damage.

A drawback of (14) is that a large number of fatigue constraints is considered. This leads to high computational effort in solving the optimization problem. Hence, the n_e fatigue constraints are replaced by a single bound, i.e.

$$ \underset{e \in \{1, 2, ..., n_{e}\}}{\max}D_{N,e}({\boldsymbol{x}}) \le \bar{D}, $$

(15)

with $\bar {D}$ as the maximum allowable damage. Since the max-function is non-differentiable, we then replace it with a smooth aggregation function (Kennedy and Hicken 2015). The left-hand side in (15) is thus approximated using the p-norm as

$$ \begin{aligned} D^{{PN}}({\boldsymbol{x}}) = \left[\sum\limits_{e=1}^{n_{e}}(D_{N,e}({\boldsymbol{x}}))^{P} \right]^{\frac{1}{P}} \end{aligned} $$

(16)

with P > 1 and D^PN as the approximated maximum damage. It holds that $D^{PN}({\boldsymbol {x}}) \to \max \limits _{e} D_{N,e}({\boldsymbol {x}})$ when $P \to \infty$, but too large values for P can cause numerical difficulties.

The mass minimization problem is

$$ \begin{aligned} (\mathbb{T}\mathbb{O}) \begin{cases} \underset{{\boldsymbol{x}}}{\min}\ \overset{n_{e}}{\underset{e = 1}{\sum}}m_{e}\rho_{e}({\boldsymbol{x}})\\ \text{s.t.}\begin{cases} D^{PN}({\boldsymbol{x}})\le\bar{D}, \\ {\boldsymbol{F}}^{\text{T}}{\boldsymbol{u}}({\boldsymbol{x}})\le\bar{C}, \\ \epsilon\le x_{e} \le 1, \quad e = 1, ..., n_{d},\\ x_{e} = \epsilon, \quad e = n_{d} + 1, ..., n_{e}. \end{cases} \end{cases} \end{aligned} $$

where m_e is the mass of each FE, n_d is the total number of elements that are considered in the design domain and $\bar {C}$ is the maximum structural compliance. The compliance constraint is included to avoid all-void solutions, see the remark in Suresh et al. (2020).

In the above optimization problem, we have implemented the domain extension approach from Clausen and Andreassen (2017), where the external boundaries are treated in the same way as the internal boundaries. In this method, the optimization variables are prescribed to the lower limit 𝜖 in the extended parts, which relate to free boundaries, while we solve the elasticity problem for the whole domain. In ($\mathbb {T}\mathbb {O}$), we assume the last set of n_e − n_d elements as the extended elements. The domain extension approach is used to get a smoother profile at the boundary and to avoid artificial structural reinforcement at the domain boundary.

We use the adjoint method to compute the derivative of the fatigue response with respect to the design variables. The predicted fatigue damage for the proposed model has history dependence so the adjoint variables are obtained by solving a (discrete) terminal value problem (c.f. the expressions in Suresh et al. 2020).

4 Numerical examples

The proposed method is implemented in the in-house FE program TRINITAS. The optimization problem ($\mathbb {T}\mathbb {O}$) is solved by using the Method of Moving Asymptotes (MMA) (Svanberg 1987), on a desktop computer with an Intel®; Core™ i7-7500U CPU @ 2.70 GHz with 24 GB of RAM.

We consider two materials, namely AISI-SAE 4340 alloy steel and 34CrMo6 forged steel. As mentioned earlier, the fatigue damage D(t) increases from D = 0 to D = 1 and with these conditions, the fatigue model parameters are calibrated against Wöhler curves following the fitting procedure described in Ottosen et al. (2008). The fitted fatigue parameters for these materials are shown in the Table 1. The AISI-SAE 4340 steel has isotropic fatigue properties (S_∥ = S_⊥ and A_∥ = A_⊥), while the 34CrMo6 forged steel has directional dependent fatigue properties, i.e. the longitudinal fatigue strength is greater than the transverse fatigue strength (S_∥ > S_⊥). In AM, when forged steel is used, the fatigue strength is sensitive to load and build (longitudinal) directions. If the load and build directions are the same, the fatigue strength is greatest, while the fatigue strength is the weakest when the directions are perpendicular to each other. For instance, this can be easily seen for a simple uniaxial case. The parameter ξ in (5) becomes ξ = 1 when the uniaxial load is applied in the direction of build direction, and as a result, the endurance stress $\tilde {S}_{0}$ is interpolated to $\tilde {S}_{0} = S_{\parallel }$ and thus have a higher fatigue strength.

Table 1 Fitted fatigue material parameters

Using the material parameters in Table 1, several examples, which are all discretized by bi-linear quadrilateral plane stress FEs having thickness 1 mm, are tested. We take the values of the penalization parameters as q = 3 and r = 0.5. The initial design variables are taken as x_e = 0.85 and the lower bound on the design variables is 𝜖 = 0.001. The problems are solved with the density filter radius as 1.5 times the element size and the exponent of the p-norm P = 8.

We notice that the p-norm function (16) overestimated the maximum value, and the adaptive scaling strategy from Le et al. (2010) and Zhang et al. (2019) may therefore be useful. However, when we compare the optimization result of the T-shaped beam (for the case b) using the p-norm (Table 3) with those from the adaptive scaling strategy for the parameter setting of Le et al. (2010) and Zhang et al. (2019), no difference in topologies were observed by the eye. Nor did we observe any such difference when comparing the topologies of the L-shaped beam (case of 90^∘ build direction) with non-proportional load (Table 4). Therefore, the results in the following are for a fixed p-norm.

The stopping criterion used for the optimization is given by

$$ \begin{array}{@{}rcl@{}} \left|\frac{f_{k+1} - f_{k}}{f_{k+1}}\right| \le \text{tol}, \end{array} $$

where f_k is the objective value at the k th iteration. The tolerance is set to tol = 0.1E − 6.

To determine the maximum structural compliance $\bar {C}$ in $(\mathbb {T}\mathbb {O})$, we perform a simple compliance minimization subjected to a mass constraint (without fatigue constraint). The optimization problem is solved for 40% of the original structural mass and the optimal compliance obtained is set as $\bar {C}$. The maximum fatigue damage $\bar {D}$ in $(\mathbb {T}\mathbb {O})$ is selected such that we try to obtain an infinite life in the structure. From Ottosen et al. (2008), Holopainen et al. (2016), and Lindström et al. (2020), it is found that the damage per cycle for getting an infinite life is around 1E − 7. Hence, for the presented examples, the maximum fatigue damage $\bar {D}$ is selected close to that value.

4.1 T-shaped beam with non-periodic proportional load history

Figure 2 shows a T-shaped beam, where L₁ = 100 mm. Owing to the symmetric profile of this geometry, it serves as a good benchmark test for solving the optimization problem with materials having isotropic or anisotropic properties. On solving the problem, we expect a symmetric profile for the optimized design when an isotropic material is used, while a non-symmetric profile is expected for an anisotropic material.

Fig. 2

Geometry of the T-shaped beam

The model is discretized by 8800 FEs, and the top edge of the beam is clamped. To account for stiffness and fatigue, two load cases are created. For the first load case, a static load Q₁ = 1200 N is applied at each one of two opposing ends of the beam, and used to evaluate the compliance. In the second load case, a Gaussian-random load, shown in Fig. 3, i.e. $\tilde {Q}_{1}(t) = Q_{1}S_{f}(t)$, with t = 0,0.05,...50, is applied for fatigue estimation. Using (12) and (13), the stresses ${\boldsymbol {\sigma }}({\boldsymbol {x}},\tilde {Q}_{1}(t))$ are first computed and then these stresses are used to estimate the fatigue damage. The grey regions in Fig. 2 indicate elements where x_e = 𝜖 for the domain extension approach. The red region below the loads contains elements that are fixed to x_e = 1 and are not subject to optimization due to unavoidable stress concentrations.

Fig. 3

Non-periodic load history

For the T-shaped beam, we solve $(\mathbb {T}\mathbb {O})$ for AISI-SAE 4340 alloy steel and 34CrMo6 forged steel, respectively. The fatigue parameters for these materials are given in Table 1. As mentioned earlier, the AISI-SAE 4340 steel has isotropic fatigue properties while 34CrMo6 forged steel has anisotropic fatigue properties. The problem is solved for various cases, shown in Table 2 along with corresponding bounds $\bar {D}$ on fatigue. In the isotropic case, the Young’s modulus is taken as 210 GPa and Poisson’s ratio as 0.3. In the transversely isotropic case, the five independent material parameters used in the material matrix (see Appendix) are E₁ = 200 GPa, E₃ = 120 GPa, G₁₃ = 78 GPa, ν₁₂ = 0.3 and ν₁₃ = 0.23. Direction 3 corresponds to the build direction defined by the unit vector m. The upper bound on compliance in the isotropic elastic cases is taken as $\bar {C} = 1.45$ Nmm, while in anisotropic elastic cases, the upper bound on compliance is $\bar {C} = 1.75$ Nmm.

Table 2 Upper bound on fatigue damage $\bar {D}$ in units of 1E − 5 for various cases; the fatigue constraint cases (b) and (e), use the fatigue parameters for AISI-SAE 4340 while the cases (c) and (f) use the fatigue parameters for 34CrMo6

Table 3 gives optimization results of the T-shaped beam for the cases from Table 2 along with the corresponding objective values. The first row provides different optimization profiles of the beam when an isotropic linear elastic model is used. In the case where $(\mathbb {T}\mathbb {O})$ is solved without fatigue constraint, the optimized design has a symmetric profile with sharp corners with high stress concentrations at the re-entrant corners of the design domain. In the cases with fatigue constraint, for both materials, the profiles obtained in the optimized designs have beam-like members with joints formed away from the re-entrant corners to reduce the high stress concentrations and thereby prolonging the structural life. With AISI-SAE 4340 steel, the material has isotropic fatigue properties (see Table 1), and hence, the design has a symmetric profile. In contrast, the optimized model with 34CrMo6 steel has a non-symmetric profile as the material has anisotropic fatigue properties. We can see that more material is distributed in the right side of the optimized profile as there is high fatigue damage distribution in the right side.

Table 3 Optimization results of the T-shaped beam for different cases showing the topology of optimized designs and then mass values. x₁ and x₃ indicate the direction that take the Young’s modulus values E₁ and E₃, respectively

The bottom row in Table 3 provides optimized designs when anisotropic material is used in the constitutive response. Similar to the top row, the optimized design obtained after solving $(\mathbb {T}\mathbb {O})$ without fatigue constraint has sharp corners at the re-entrant corners of the design domain, but the profile is non-symmetric due to the anisotropic constitutive response. When a fatigue constraint is used, we obtain, for both materials, optimized designs with profiles having beam-like members with joints formed away from the re-entrant corner to reduce the high stress concentrations. Unlike the previous case, the optimized design obtained with AISI-SAE 4340 steel has a non-symmetric profile. This is due to the anisotropic constitutive response even if the material has isotropic fatigue properties. With 34CrMo6 steel, the optimized model has a non-symmetric profile due to the combination of anisotropic constitutive response and anisotropic fatigue properties of the material.

The computational effort for solving the above examples with fatigue constraint is quite high. The times required to get the optimized designs in Table 3 are around 1.5 h for (b) (converged after 210 iterations); 4 h for (c) (converged after 680 iterations); 4 h for (e) (converged after 720 iterations); and 3.5 h for (f) (converged after 542 iterations). Most time is spent in sensitivity analysis as it depends not only on total number of elements but also on the number of time steps.

4.2 L-shaped beam with non-proportional load history

The L-shaped beam is shown in Fig. 4, with L₂ = 100 mm. The design domain is discretized by 6400 FEs and the top edge is clamped. Like the previous example, two load cases are created. For the first load case we apply static loads Q₂ = 700 N, Q₃ = 1000 N and Q₄ = 500 N for compliance evaluation. In the second load case, we apply an arbitrary load history, shown in Fig. 5, i.e. $\tilde {Q}_{2}(t) = Q_{2}S_{1}(t)$, $\tilde {Q}_{3}(t) = Q_{3}S_{2}(t)$ and $\tilde {Q}_{4}(t) = Q_{4}S_{3}(t)$, with t = 0,0.1,...,25 for fatigue estimation. Similarly to the previous examples, the grey region contains elements where x_e = 𝜖 and the red region contains the elements that are excluded from optimization due to unavoidable stress concentrations. As a remark we note that a common approach when optimizing the L-shaped beam is to use a square mesh but fixing the design variables in the upper right corner to their lower bound to simulate an L-shaped domain. This effectively amounts to use boundary extension on the upper right parts.

Fig. 4

Geometry of the model

Fig. 5

Arbitrary load histories applied to Q₂, Q₃ and Q₄, respectively

With anisotropic constitutive response, the fatigue constrained problem $(\mathbb {T}\mathbb {O})$ is solved for the anisotropic, 34CrMo6 forged steel. In Kumara et al. (2018), anisotropic elastic properties in alloy 718 is modelled. It was observed that the value of Young’s modulus in build direction is about 70–75% of the Young’s modulus value normal to the build direction. Hence, the five independent material parameters used in the material matrix of 34CrMo6 steel are E₁ = 200 GPa, E₃ = 150 GPa, G₁₃ = 78 GPa, ν₁₂ = 0.3 and ν₁₃ = 0.23, with direction 3 being the build direction. The problem is solved for the two build directions 0^∘ and 90^∘, and the maximum fatigue damage for each build direction is $\bar {D} = 0.8\text {E}-5$ and $\bar {D} = 0.2\text {E}-5$, respectively. The upper bounds on compliance for the respective directions are $\bar {C} = 3.5$ Nmm and $\bar {C} = 4.0$ Nmm.

Table 4 provides optimization results of the L-shaped beam along with the corresponding objective values. The first row of the table provides the topology of the optimized model for each direction, while the second row gives the fatigue damage distribution within the structure. We notice a considerable change in the design for each build direction. The static loads Q₂ and Q₃ induce a high bending moment in the region close to the clamped edge. Therefore, on optimization, fatigue damage starts to dominate at the vertical boundaries near the clamped edge due to high bending stresses. For the case of 0^∘ build direction, the surface normal to the unit vector m, i.e. the vertical surface close to the clamped edge, is highly sensitive to the fatigue damage, as seen in the fatigue damage plot. Hence, more material is distributed near this surface. In the case of 90^∘ build direction, horizontal surfaces are sensitive to fatigue damage. Hence, we notice a more uniform fatigue damage distribution and also a more uniform distribution of material within the structure compared with the damage in the optimized structure obtained from 0^∘ build direction.

Table 4 Optimization results of the L-shaped beam with different build directions. x₁ and x₃ indicate the direction that take the Young’s modulus values E₁ and E₃, respectively. The first row gives the topology of the optimized models and the second row provides the fatigue, D_N,e in units of 1E − 5 within the structure

The evolution of the objective function and the fatigue constraint for the case of 90^∘ build direction is shown in Fig. 6. We notice that there is some fluctuation of the plot in early iterations but not close to convergence. The fluctuations are related to updates of the asymptotes in the MMA method. The computation time required to solve the above cases is around 4.5 h (converged after 543 iterations) and 6 h (converged after 760 iterations), respectively.

Fig. 6

Convergence plots for the case of 90^∘ build direction: Left: Objective function; Right: Fatigue constraint

The optimized design of the L-shaped beam is not very sensitive to the particulars of the statistical realization of the load. To demonstrate this, we solve ($\mathbb {T}\mathbb {O}$) by applying the load histories shown in Fig. 5, but with different amplitude values, i.e. $\tilde {Q}_{2}(t) = Q_{2}S_{2}(t)$, $\tilde {Q}_{3}(t) = Q_{3}S_{3}(t)$ and $\tilde {Q}_{4}(t) = Q_{4}S_{1}(t)$, with t = 0,0.1,...,25 for fatigue estimation. The optimized result obtained in the case 90^∘ build direction has a similar topology to the previous case shown in Table 4. This implies that when multiple statistical realizations of the load is applied to the L-shaped geometry, we can expect a similar optimized design profiles as in Table 4.

Since we employ a design variable filter in the above numerical examples, there exist grey transition regions in the design domain where the design variables have intermediate values between 0 and 1. These regions may be eliminated or diminished by using, e.g. the Heaviside projection filter by Guest et al. (2004) or by the two-step strategy in Holmberg et al. (2015). However, for simplicity this is not used in the examples presented in the paper.

4.3 3D L-shaped beam with non-periodic load history

For the final example, we take a 3D L-shaped beam, where the dimensions of the geometry are shown in Fig. 7, with L₃ = 100 mm. The design domain is discretized by 6720 eight-noded trilinear brick elements. Two load cased are created. For the first load case, we apply static load Q₅ = 1E6 N/mm and for the second load case, a non-periodic load, shown in Fig. 3, i.e. $\tilde {Q}_{5}(t) = Q_{5}S_{f}(t)$, with t = 0,0.05,...50, is applied for fatigue estimation. The grey regions in Fig. 7 indicate elements where x_e = 𝜖 for the domain extension approach.

Fig. 7

Geometry of the 3D L-shaped beam

The fatigue constrained problem $(\mathbb {T}\mathbb {O})$ is solved for the pure isotropic case with AISI-SAE 4340 steel. In the isotropic elastic response, the Young’s modulus is taken as 210 GPa and Poisson’s ratio as 0.3. The maximum fatigue damage in $(\mathbb {T}\mathbb {O})$ is $\bar {D} = 0.1\text {E}-5$, while the upper bound on compliance is $\bar {C} = 10$ Nmm.

Figure 8 provides the optimized result of the 3D L-shaped beam. Here we present an iso-surface based on ρ = 0.5 of the optimized design. The final mass of the optimized design after 862 iterations is 0.19 kg (≈ 47% of the original mass) and computation time required to solve this problem is around 9 h.

Fig. 8

Optimized model of the 3D L-shaped beam

5 Conclusion

Using the continuous-time HC fatigue model, the gradient-based TO formulation with HC fatigue constraint (Suresh et al. 2020) was extended to handle transversely isotropic material properties. The presented approach is capable of handling arbitrary load histories, including non-proportional load histories.

Fatigue sensitivities were derived by the adjoint method. The domain extension approach was implemented to get a smoother profile at the boundary. The proposed method is capable of handling 2D and 3D problems. However, the computational effort is quite high, mainly due to sensitivity analysis. Hence, an interesting extension would be to use acceleration techniques to reduce the computational time.

Several numerical examples were given. First, a T-shaped beam subject to a non-periodic load history which serves as a comparison test between isotropic and anisotropic materials, in this case AISI-SAE 4340 alloy steel and 34CrMo6 forged steel. For the isotropic material with isotropic fatigue properties, the optimized design has symmetric profiles whereas for the anisotropic material (even with isotropic fatigue properties), the designs have non-symmetric profiles. The second numerical example was an L-shaped beam subject to non-proportional loads tested with two build orientations, 0^∘ and 90^∘. For the 0^∘ build direction, the material within the structure is heavily distributed in regions near the vertical surfaces, while the 90^∘ build direction yields more uniform distribution of material within the structure. Finally, we also presented a 3D L-shaped beam with non-periodic history for a pure isotropic case, thus demonstrating the capability to handle 3D problems.

Appendix: Transversely isotropic material

For a transversely isotropic material, the material matrix E has five independent parameters. In the longitudinal direction, i.e. direction 3, the Young’s modulus and Poisson’s ratio are E₃ and ν₁₃, while in the transverse plane, i.e. directions 1 and 2, the material has a symmetry plane with E₁ = E₂ and ν₁₂ as Young’s modulus and Poisson’s ratio, respectively.

The compliance matrix, which is the inverse of the material matrix, is written as

$$ {\boldsymbol{E}}^{-1} = \left[ \begin{array}{c c c c c c} \frac{1}{E_{1}} & \frac{-\nu_{12}}{E_{1}} & \frac{-\nu_{13}}{E_{1}} & 0 & 0 & 0 \\ \frac{-\nu_{12}}{E_{1}} & \frac{1}{E_{1}} & \frac{-\nu_{13}}{E_{1}} & 0 & 0 & 0 \\ \frac{-\nu_{13}}{E_{1}} & \frac{-\nu_{13}}{E_{1}} & \frac{1}{E_{3}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{13}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{13}} \end{array}\right], $$

where the shear modulus in the symmetry plane is $G_{12} = \frac {E_{1}}{2(1 + \nu _{12})}$.

From the compliance matrix, the five independent parameters are E₁, E₃, G₁₃, ν₁₂ and ν₁₃. In the case of isotropic materials, we have E₁ = E₃, ν₁₂ = ν₁₃ and G₁₂ = G₁₃.