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

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.


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 AMfabricated 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  which use cycle-counting algorithms. Other examples of structural optimization with fatigue constraints are found in , 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 nonproportional loads and is implemented in TO problems. Our recent contribution  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  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 pnorm (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.

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 continuoustime 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. β > 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 aṡ where C > 0, K > 0 and L > 0 are material parameters.
The deviatoric stress is s = σ − 1 3 tr(σ )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.

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.
where the transverse component σ ⊥ is calculated by the use of a projection tensor P defined as P = I − M. This gives To formulate the endurance function we introduce four tensor invariants, Following Holopainen et al. (2016), the endurance function is written as where A > 0 and A ⊥ > 0 are material parameters related to the corresponding physical directions. It can be shown that tr(σ ) = I 4 1 and tr(σ ⊥ ) = I 1 − I 4 . The effective stress σ is obtained bȳ The scalar function ξ is the stress ratio, i.e.
Taking the time derivative of (4) and with help of (5), we geṫ Substituting (1) forα, and rearranging, we have When H (β) > 0 and since S ⊥ + H (β)H (β)Cσ > 0, we can write the above rate function aṡ The isotropic HCF model in Ottosen et al. (2008) emerges as a special case when A ⊥ = A = A > 0 and S ⊥ = S = S 0 > 0. Then the endurance function (4) and its rate (6) become With the endurance function and its rate defined in (4) and (6), the fatigue damage is estimated by integrating the ODEs (1) and (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 The discrete versions of (1) and (2) are written as The increment of the endurance function in (6) is approximated as with the stress increment as σ i = σ i − σ i−1 , and the effective stress asσ i =σ (σ i , α i ). The invariants are computed as 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.

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-andwhite designs. A design variable filter (Bruns and Tortorelli 2001) is applied to avoid mesh dependency and checkerboard patterns arising for commonly used low-order FEs. This filter gives a new set of variables ρ = ρ(x) = [ρ 1 (x), ρ 2 (x), ..., ρ n e (x)] T called physical variables since they define the structural stiffness and the mass.
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 historyF (t) is used for estimating the fatigue damage.
The equilibrium equation for the static load case is For the fatigue load case, we assume quasi-static equilibrium, and consider a sequence of static problems, In the equilibrium equations, K(ρ(x)) is the global stiffness matrix, and u andũ(t) are the displacement vectors for the corresponding load cases. At each time step t i = i t, we have whereF i =F (t i ), and the solution for a given design x and time step i is denoted byũ i (x).
To promote black-and-white designs, the SIMP approach (Bendsoe and Sigmund 2004;Christensen and Klarbring 2009) is used, i.e.

K(ρ(x))
where K e are elemental stiffness matrices and q > 1 is a penalization parameter. For a givenũ i (x), the stress vector at a stress evaluation point is calculated aŝ 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 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 T u(x). As for the fatigue load case, for a given load historyF 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 withD 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.
withD as the maximum allowable damage. Since the maxfunction 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 with P > 1 and D P N as the approximated maximum damage. It holds that D P N (x) → max e D N,e (x) when P → ∞, but too large values for P can cause numerical difficulties.
The mass minimization problem is 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 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 (TO), 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).

Numerical examples
The proposed method is implemented in the in-house FE program TRINITAS. The optimization problem (TO) 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 stressS 0 is interpolated toS 0 = S and thus have a higher fatigue strength.
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 where f k is the objective value at the kth iteration. The tolerance is set to tol = 0.1E−6.
To determine the maximum structural complianceC in (TO), 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 asC. The maximum fatigue damageD in (TO) is selected such that we try to obtain an infinite life in the structure.  Figure 2 shows a T-shaped beam, where L 1 = 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.

T-shaped beam with non-periodic proportional load history
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 1 = 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.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 σ (x,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.
For the T-shaped beam, we solve (TO) 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 boundsD 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 1 = 200 GPa, E 3 = 120 GPa, G 13 = 78 GPa, ν 12 = 0.3 and ν 13 = 0.23. Direction 3 corresponds to the build   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 (TO) 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. 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 (TO) 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.

L-shaped beam with non-proportional load history
The L-shaped beam is shown in Fig. 4, with L 2 = 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 2 = 700 N, Q 3 = 1000 N and Q 4 = 500 N for compliance evaluation. In the second load case, we apply an arbitrary load history, shown in Fig. 5, i.e. Q 2 (t) = Q 2 S 1 (t),Q 3 (t) = Q 3 S 2 (t) andQ 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.
With anisotropic constitutive response, the fatigue constrained problem (TO) 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 Table 3 Optimization results of the T-shaped beam for different cases showing the topology of optimized designs and then mass values. x 1 and x 3 indicate the direction that take the Young's modulus values E 1 and E 3 , respectively direction. Hence, the five independent material parameters used in the material matrix of 34CrMo6 steel are E 1 = 200 GPa, E 3 = 150 GPa, G 13 = 78 GPa, ν 12 = 0.3 and ν 13 = 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 isD = 0.8E−5 andD = 0.2E−5, respectively. The upper bounds on compliance for the respective directions areC = 3.5 Nmm andC = 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 2 and Q 3 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.
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.
The optimized design of the L-shaped beam is not very sensitive to the particulars of the statistical realization of  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 Table 4 Optimization results of the L-shaped beam with different build directions. x 1 and x 3 indicate the direction that take the Young's modulus values E 1 and E 3 , 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  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.

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 3 = 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 5 = 1E6 N/mm and for the second load case, a The fatigue constrained problem (TO) 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 (TO) isD = 0.1E − 5, while the upper bound on compliance isC = 10 Nmm. Figure 8 provides the optimized result of the 3D Lshaped 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. Using the continuous-time HC fatigue model, the gradientbased TO formulation with HC fatigue constraint ) 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.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.

Replication of results
To reproduce the above optimization results, the fatigue model in Section 2 is first implemented. The optimization problem (TO) is solved by the MMA method from Svanberg (1987). The fatigue sensitivities are derived by the adjoint method and the procedure is given in Suresh et al. (2020). Relevant parameters are given in Section 4 to run the numerical examples.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.

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 3 and ν 13 , while in the transverse plane, i.e. directions 1 and 2, the material has a symmetry plane with E 1 = E 2 and ν 12 as Young's modulus and Poisson's ratio, respectively.
The compliance matrix, which is the inverse of the material matrix, is written as where the shear modulus in the symmetry plane is G 12 = E 1 2(1 + ν 12 ) .
From the compliance matrix, the five independent parameters are E 1 , E 3 , G 13 , ν 12 and ν 13 . In the case of isotropic materials, we have E 1 = E 3 , ν 12 = ν 13 and G 12 = G 13 .