Topology optimization accounting for surface layer effects

Metal AM (additive manufacturing) components are generally inhomogeneous and have different microstructure in the bulk compared with (contour) regions near the surface. This, as well as rough as-built surfaces, affects mechanical properties. In this paper, we develop a topology optimization method that considers such inhomogeneities. The method is a direct extension of standard density-based methods using linear filtering for regularization, and a second filtering of the design variables is used to identify a surface layer, the thickness of which is given by the filter radius. Domain extension is used in order to properly identify such layers at the boundary of the design domain. The method is generally applicable but is demonstrated for stiffness optimization. Both two- and three-dimensional problems are treated. A general property of the method is that the topological complexity is reduced, i.e. the optimized designs get fewer and thicker structural members as the width of the surface layer is increased.


Introduction
Topology optimization (TO) and additive manufacturing (AM) is an unusually natural and powerful match of design and manufacturing methods, sharing the property of very large freedom in geometrical form. This fact has prompted new developments in TO, tailored to the needs of AM, such as overhang constraints and anisotropic material properties. However, the important fact that mechanical properties of metal AM structures frequently show a pronounced member size dependency (Algardh et al. 2016), originating from surface or boundary layer effects, has not been modelled and incorporated in TO methods. The present paper develops a TO method that identifies boundary regions and assigns unique properties, different from bulk properties, to these regions. Kahlin et al. (2017) studied the mechanical properties of AM titanium Ti-6Al-4V test specimens built by both Responsible Editor: Ji-Hong Zhu Anders Klarbring anders.klarbring@liu.se 1 Division of Solid Mechanics Department of Management and Engineering Institute of Technology, Linköping University, 581 83, Linköping, Sweden electron beam melting and laser sintering. In particular, they compared fatigue properties of specimens that have rough as-built surfaces with specimens where surfaces have been machined and polished. A distinct trend in these results is that the fatigue limit is increased with improved surface quality. Moreover, even when surface roughness effects are essentially removed by surface treatments, there is still a difference in material properties in the bulk compared with the boundary region or surface layer. Balachandramuthi et al. (2019) studied the microstructure of AM alloy 718, investigating specimens made by both electron beam melting and selective laser sintering, and found that the hatch (bulk) region has pronounced directional properties, resulting in anisotropic material behaviour, while the contour (surface layer) region shows much less directionality in the microstructure, indicating an isotropic material behaviour. Anisotropic material behaviour of metal AM components is quantified in Kumare et al. (2018). In order to take into account both contour region effects and surface roughness effects, we develop a particular TO method as indicated above.
TO methods that identify boundary regions have previously been developed for application to coated structures, where the coating material typically has a higher stiffness than the bulk material. This represents the inverse of the present physical motivation, where, e.g. a rough asbuilt surface implies a boundary region of less stiffness then the bulk. TO methods for coated structures have been developed based on both density methods (Clausen et al. 2015;Yoon and Yi 2019;Wu et al. 2017) and level-set methods (Wang and Kang 2018). These density methods use a sequential application of filtering and projections, which is in contrast to the more direct approach in this paper, where two separate densities calculated by standard linear filtering are used, as shown in Fig. 1.

Filtering and material interpolation
Motivated by the experimental studies discussed in the introduction, we develop a TO method that identifies surface layers of the optimized structure and attributes different material properties to bulk and surface. Three density-like field variables are used: • The optimization variable ξ . • The shape density ρ, obtained from ξ by standard linear filtering. • The surface layer identifier ρ S , also obtained from ξ by standard linear filtering, but for a large filter radius than for ρ.
The first two of these variables are used in a way that is standard in modern TO Bendsøe and Sigmund (2002) and Christensen and Klarbring (2009), while the surface layer identifier ρ S is novel and used to calculate a zero-one-valued function where H is the Heaviside step function, ∇ is the gradient operator and | · | denotes the Euclidean norm. Figure 1 illustrates the relation between these functions in the ideal situation when the optimization variable ξ is zero-onevalued. The function S-which at an optimum of a TO problem is expected to take the value zero in the bulk of the structure and the value one close to the boundary-can be used to interpolate essentially any material parameter of interest. In this paper, we study two cases: Firstly, a simple model where different but constant Young's moduli in the bulk and in the surface layer are considered. These are denoted E B and E S , respectively, and the interpolation of Young's modulus, using SIMP (Bendsøe and Sigmund 2002;Christensen and Klarbring 2009) for the shape density, is then where q ≥ 1 is the SIMP exponent. Secondly, we interpolate between transversely isotropic material in the bulk, represented by a matrix E B , using Voigt notation, and isotropic material in the surface layer, where the matrix is denoted E S . Expression (2) is then replaced by Equation (2) can be included as a special case of (3), and in the following we use the latter. The precise form of the material matrix E for both isotropic and transversely isotropic materials is given in Appendix 1.  Fig. 1 Two linear filter operations, with different filter radii R and R S , R ≤ R S , are applied to a zero-one-valued optimization variable ξ , producing ρ and ρ S . The latter is used to calculate the zero-one-valued function S in (1)

Linear filtering and gradient
This subsection introduces the linear filter operator (Borrvall 2001;Bourdin 2001) and gives the details of an approximate calculation of the gradient ∇ρ S to be used in a finite element method.
The linear filter operator takes the field ξ , defined on the design domain , and delivers the field ρ S , defined on the same domain, by the convolution The integral kernel ψ is given by where R S is a filter radius. This definition implies ψ(x, y) dy = 1, for all x ∈ , The field ρ is calculated from ξ in the same way, but using a filter radius R ≤ R S . We like to calculate the Euclidean norm of the gradient of ρ S (x), and do so by first calculating partial derivatives. Since the derivative ∂ψ/∂x j is defined almost everywhere, we have ∂ρ S ∂x j (x) = ∂ψ ∂x j (x, y)ξ(y) dy.
For x and y such that ∂ψ/∂x j exists and is non-zero, i.e. when x = y and |x − y| < R S , we have where P j (y − x) is the projection of the unit vector (y − x)/|y − x| on the j -axis, i.e. cos or sin of an angle, and x = {y ∈ : |x − y| < R S }. Note that for points x a distance R S away from the boundary of , F (x) is independent of x and equal to (π R d+1 S )/3, where d is 2 or 3 depending on the spatial dimension of the problem.

Discrete filter derivative
Equation (5) is discretized in order to be used in a finite element context. Taking ξ as element-wise constant, having a value ξ e for finite element e, we get, where e is the domain of element e. Approximating further by taking one integration point per element, we get where y e C is the centroid (integration point) of element e and the sum needs to be taken only over those elements where ∂ψ/∂x j is non-zero. This equation is evaluated at x = x i C , the centroid of element i, and the other functions in the expression of ∂ψ/∂x j are approximated as Having calculated the partial derivatives of ρ S , we have The discretized zero-one measure used in the optimization problem, defined in the next section, is then where the smooth functioñ is used to approximate the Heaviside step function for x ≥ 0 by choosing A large enough. For illustration, approximate expressions for |∇ρ S (x)| 2 in case of a uniform 2D mesh are given in Appendix 2 for two different ranges of the filter radius. Note, however, that the numerical implementation uses the general formulas (6) and (7), valid for arbitrary meshes.

Discrete optimization problem
Given a standard displacement-based finite element discretization of linear elasticity, we consider a stiffness optimization problem as follows: Here, n is the number of finite elements covering , V is the total available structural volume and a i is the volume of element i. The vector ξ T = (ξ 1 , . . . , ξ n ) is the discretized optimization variable, where we use the standard setting of one such variable for each finite element. The corresponding discrete shape density ρ is a function of ξ through the linear filter defined in the previous section. Following Clausen and Andreassen (2017), we use the domain extension approach in order to treat external boundaries in the same way as internal boundaries. This is crucial for identification of boundary regions by means of the filtered density ρ S . If is considered the extended domain, the domain extension method means that optimization variables are prescribed to the lower box limit in extended parts (shown in Figs. 2, 4 and 5), while the elasticity problem is solved for the whole domain . Without loss of generality, we assume that the last n − m elements of ξ are fixed to the lower bound .
For a given external load vector F , u = u(ξ ) is the solution of the equilibrium state problem where using (3) the stiffness matrix can be written as where K S i and K B i are (global) element stiffness matrices of element i for surface and bulk material, respectively, S i (ξ ) was defined in Section 2.1 and the function ρ i (ξ ) is given by a discretized version of the linear filter operator (Borrvall 2001), defined as where the filter weights are The small positive number in (SO), together with suitable support conditions, guarantees that K(ξ ) is nonsingular for all admissible designs, giving a unique solution u = u(ξ ). Problem (SO) is solved by the method of moving asymptotes (MMA) (Svanberg 1987) using expressions for the sensitivity of the objective function given in Appendix 3. The computational cost is dominated by solving (9), similar to standard stiffness based TO. The stopping criterion for the optimization is given by where f k is the objective value at the k th iteration and tol is a given value.
Remark 1 As mentioned in the introduction, density-based TO methods that identify boundary regions have previously been developed for application to coating material (Clausen et al. 2015;Yoon and Yi 2019;Wu et al. 2017). The comparison with the present approach is, however, not direct since the aim is a different application and since the methods differ in several important ways. The present approach uses two separate applications of linear filtering of the design variable followed by one projection (the Heaviside step), while Clausen et al. (2015) use a sequence of filterings and projections, where continuation strategies are needed for the projection parameters. Moreover, to enable calculation of gradients, the filters in Clausen et al. (2015) are of PDE-type using nodal-based densities. The present work shows that gradients can be effectively calculated for classical elementwise densities by using the basic analytical expression (4) for the linear filter. The actual formulation of the optimization problem also differs between Clausen et al. (2015) and the present work. The volume constraint of problem (SO) is such that the "cost" is the same for bulk and surface layer material. In Clausen et al. (2015), the corresponding constraint is formulated in terms of mass and by giving different physical densities for bulk and surface layer material, the "cost" will be different for the two types of material. Moreover, they also use a different interpolation from (2), which in the notation of this paper is It is not clear at the present stage of development whether the present method, intended for AM applications, could be used also for the application to coated structures and vice versa.

Numerical results
The proposed method is implemented in the in-house finite element program TRINITAS (2000). Results are presented for both the isotropic case, when (2) is used for material interpolation, and for the case of transversely isotropic bulk material, using (3). All calculations are performed on a desktop computer with an Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz with 24GB of RAM. The equilibrium state problem (9) is solved with a sparse direct solver from Intel MKL.
For the case of isotropic material properties, we use a bulk material having Young's modulus 210 GPa and Poisson's ratio 0.3. This Young's modulus is then lowered in the surface region to represent, e.g. the influence of the rough as-built surface. Presently, no experimentally determined value for this surface region Young's modulus is known, but the qualitative influence of such an inhomogeneity on optimal topologies is certainly of interest and in the examples we use E B = 210 GPa and E S = 105 GPa. For the case of transversely isotropic bulk material, we use (3) to interpolate between transversely isotropic material in the bulk, represented by E B , with the independent material constants defined in Appendix 1 set to E 1 = 210 The value of the SIMP parameter is q = 3 and the initial design variables are taken as ξ e = 0.5. The lower bound on the design variable is ε = 0.001. The maximum bound on the structural volume is V = 0.4V 0 in the 2D problems and V = 0.3V 0 in the 3D problem, with V 0 as the initial volume of the structure. The parameter A in the Heaviside step function approximation should be scaled with respect to finite element size h and filter radius R S . Since the approximate expression for |∇ρ S (x)| 2 is proportional to h 2d /R 2+2d S , see Appendix 3 for the 2D case d = 2, we use A = γ (R 2+2d S /h 2d ) where γ = 8. The stopping tolerance in (10) is set to tol = 10 −9 .

2D cantilever beam
As the first example, we consider a 2D cantilever beam, shown in Fig. 2. The dimensions of the geometry are shown in Fig. 2, where L 1 = 200 mm. The design domain is discretized by 200 × 120 bi-linear square plane stress elements with thickness of 1 mm and side length h = 1 mm. The beam is fixed at the left end, while a static load F 1 = 150 N, distributed over 4 mm, is applied at the other end. The regions in grey colour indicate elements used for domain extension.
The standard filter radius used to obtain the shape density is taken as R = 1.2 mm, while three different values are used for the second filter radius identifying the surface layer, namely R S = 1.2, 2.3 and 3.2 mm. Table 1 shows results for the cantilever beam in the case of isotropic material. The first row of the table shows the topology of the optimized model, while the second row shows the material distribution within the structure, red colour indicating bulk material with E B = 210 GPa and green colour surface material with E S = 105 GPa. Different filter radii R S , implying different surface layer thicknesses, are compared. What is evident is that a thicker layer promotes a less complicated topology, having thicker The first row shows the topology of the optimized model and the second row shows the material distribution, where E B is red and E S is green. Length (L S ) of the boundary (perimeter) and optimal compliances are shown in the bottom lines Length (L S ) of the boundary (perimeter) and optimal compliances are shown in the bottom lines structural members. The reason for this is that the total length of the boundary is a direct measure of the proportion of the available material that has low stiffness. To quantify this, we calculate the length of the boundary of the optimized geometries in Table 1, by summing the areas of all finite elements belonging to the surface layer and dividing by the layer thickness (that is the filter radius R S ). These lengths, denoted L S , are given in Table 1. For comparison, to the far left in Table 1, the optimized topology for the standard setting without the surface layer influence is shown. Table 2 shows optimized results for the cantilever beam in case of transversely isotropic material in the bulk and isotropic material in the surface layer. As for the previous case, the first row shows the topology of the structure and the second row shows the material distribution. The red colour corresponds to transversely isotropic material which is distributed in the bulk, while the green colour indicates isotropic material properties which is distributed in the surface layer. A similar trend in geometry change as the second filter radius is increased is observed in both Tables 1 and 2, i.e. the structure changes to a less complicated topology. For both cases, a gradual increase in compliance value is noticed when R S is increased from zero. Figure 3 shows the evolution of the objective function for both cases when R S = 3.2 mm. The increase in the compliance value noticed in the initial iterations is a consequence of using an infeasible point (ξ e = 0.5) as initial design: the  The first row shows the topology of the optimized model and the second row shows the material distribution with E B in red and E S in green. Length (L S ) of the boundary (perimeter) and optimal compliances are shown in the bottom lines objective function increases until the volume constraint is satisfied. The slight oscillations noticed in the plots are due to updates of the asymptotes in the MMA method, and is not related to the particular problem formulation.

2D MBB beam
The next example is an MBB beam where the dimensions of the geometry are shown in Fig. 4, where L 2 = 100 mm. Exploiting symmetry and again using an element size h = 1 mm, the right half of the beam is discretized with 300×100 bi-linear square quadrilateral plane stress elements with a thickness of 1 mm. A static load F 2 = 150 N is applied. This force, as well as the support in the lower right corner, is distributed over 4 mm. Like in the previous example, grey colour indicates parts used in the domain extension approach.
The filter radius R = 1.35 mm is used for the shape density, while three different values are compared for the second filter radius: R S = 1.8, 2.5 and 3.5 mm. Like in the previous example, the problem is solved for both the isotropic and the transversely isotropic cases.
Results for the isotropic case are shown in Table 3. The first row shows the topology of the optimized model, while the second row shows the material distribution, with red colour indicating regions with Young's modulus E B = 210 GPa and green colour showing the surface layer with Young's modulus E S = 105 GPa. Table 4 shows results when transversely isotropic material properties are used in the bulk and isotropic material properties in the surface layer. The topology of the optimized design and the corresponding material distribution within the structure are shown in the first and second row, respectively. In all examples, it is seen that when the  Length (L S ) of the boundary (perimeter) and optimal compliances are shown in the bottom lines second filter radius is increased, the structures are changed towards less complicated topology. The length of perimeters is calculated and given in the tables. We also note that in all examples there is an increase in the optimal compliance value when the second filter radius is increased.

3D cantilever beam
For the final example, we study a 3D cantilever beam, where the dimensions of the geometry are shown in Fig. 5 with L 3 = 100 mm. The design domain, represented in orange colour, is discretized by 100 × 60 × 30 eight-noded trilinear brick elements, resulting in an element size h = 1 mm. The surface of the beam is fixed at the left end, while a static line load F 3 = 150 N/mm is applied at the other end. The regions in grey colour, shown in Fig. 5, indicate the elements that are used in the domain extension approach. The filter radius R is set to 1.2 mm and for brevity we treat only the isotropic case of (2).  Figure 6 shows the optimized design for R S = 0, i.e. no surface layer is present. The iso-surface of the optimized design, based on ρ = 0.5, is shown to the left in Fig. 6, while sliced surfaces of this design are shown to the right. The slicing is done in orthogonal directions: an XY-symmetry plane and different YZ-planes at regular intervals of 0.2L 3 . It is evident that the structure is hollow, containing an inner cavity. Figure 7 shows an optimized design for R S = 1.5 mm, with isotropic material in both bulk and surface layer. The iso-surface of the optimized topology is shown in Fig. 7. Both the geometry and the material distribution are shown using slicing. Red colour in the material distribution plot indicates bulk material having the modulus value E B = 210 GPa, while green colour indicates the surface layer with the modulus value E S = 105 GPa. We notice a change in profile Compliance is 0.25 Nmm when compared with the optimized profile obtained with R S = 0, and the inner cavity has a smaller projected area in the XY-plane. There is a change in optimal compliance value from 0.22 Nmm in Fig. 6 to 0.25 Nmm in Fig. 7, due to presence of material with low stiffness in the surface layer.

Concluding remarks
We have shown that linear filtering followed by a smoothed Heaviside step can be used to construct a function S that identifies a surface layer in TO optimized structures. The method is based on calculating the gradient of the filtered design field, using its analytical expression as a convolution integral. After a simple discretization involving elementwise constant fields and one integration point per finite element, a value S i = S i (ξ ) for each element i indicates bulk (S i ≈ 0) or boundary (S i ≈ 1). The computational cost is only minutely larger than that for standard stiffness optimization involving homogenous material, where the solution of the state problem dominates the computational time. As a general conclusion concerning the influence of a surface layer on optimized structures, it is evident that the larger the surface layer, the more topologically simple the structure becomes. The simple intuitive explanation is that for a topologically complicated structure, consisting of many structural members, the total length of the boundary is larger, resulting in more low stiffness material spent in the surface layer, working essentially like a perimeter constraint (Borrvall 2001). It can also be mentioned that similar trends have been observed in robust topology optimization considering uncertain geometrical boundary imperfections (Zhang and Kang 2017).
A goal of the presentation of the method has been to stay as close as possible to a standard setting using linear filters and element-wise densities, and, therefore, other filterrelated issues such as length-scale control and projections to achieve black-and-white designs have not been discussed. However, since the identification of the boundary layer is uncoupled, or done in parallel, to the identification of the geometrical domain, such extensions of the method are rather direct.

Replication of results
To reproduce the above optimized results, a standard finite element implementation for solving (9) is needed. The optimization problems (SO) are then solved using the MMA method from Svanberg (1987), using sensitivities derived in Appendix 3. Relevant parameters are given in Section 4.
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://creativecommons. org/licenses/by/4.0/.

Appendix 1. Transversely isotropic material
The material matrix in Voigt notation E contains 5 independent material constants for a transversely isotropic material. This is most clearly visible in the compliance matrix, which is the inverse of E. Given that direction 3 is the build direction and, thus, directions 1 and 2 define the plane of isotropy, it reads where symmetry requires ν 13 /E 1 = ν 31 /E 3 and the shear modulus in the plane of isotropy is G 12 = E 1 2(1 + ν 12 ) .
For a slightly larger filter radius, √ 2h < R S < 2h, we have where the first square parenthesis represents the derivative in the vertical (row) direction and the second square parenthesis is a similar expression representing derivative in the horizontal (column) direction.

Appendix 3. Sensitivity analysis
Since we are using the gradient-based method MMA (Svanberg 1987) for solution of (SO), we need the sensitivity or derivative of the objective function f (ξ ). The following expression is valid (Christensen and Klarbring 2009): in whichH is the derivative of the smooth Heaviside function approximation, and u denotes the solution to (9).