A nonlocal approach to graded surface modeling in topology optimization

Article Content

1 Introduction

This paper introduces a novel approach to topology optimization, focusing on incorporating surface grading within the optimized structures. Surface grading, a distinctive characteristic of additively manufactured structures, is intricately linked to both the manufacturing process and the mechanical properties inherent in the base material. Notably, in powder-bed additive manufacturing techniques such as electron beam melting (EBM) or selective laser sintering (SLS), various material properties, including elastic modulus, Poisson’s ratio, ultimate tensile strength, and porosity, exhibit non-uniformity across the thickness of the fabricated features. This non-uniform distribution leads to thickness-dependent effective material characteristics, exerting a pronounced influence on the overall properties of the structure. Such effects are particularly significant in thin-walled structures, where the thickness of the walls is comparable to the size of the laser beam spot used in the additive manufacturing process.

Multiple studies have investigated the thickness dependency of the material properties in additive manufacturing. Algardh et al. (2016) studied the thickness dependency of microstructures in thin-walled titanium parts manufactured by EBM, revealing a significant impact on mechanical properties attributed to high surface roughness. Tasch et al. (2018) determined that structures possessing a thickness less than 1 mm exhibited notable reductions in stiffness, ultimate tensile strength, and elongation at break. Sindinger et al. (2020) demonstrated evidence of thickness dependency of the material properties of laser-sintered polyamide 12 (PA12) specimens, with the degree of variation in mechanical properties changing as wall thickness decreases. Subsequently, the authors developed finite element models for the thickness-dependent Young’s modulus and Poisson’s ratio of shell structures, calibrated from experiments with laser-sintered short-fiber-reinforced PA12 material (Sindinger et al. 2021a, b, c). The thickness dependency was modeled using polynomial functions of thickness for Young’s modulus and Poisson’s ratio. On the other hand, Jaksch et al. (2022) introduced a nonlocal integral model to incorporate surface grading of Young’s modulus, simulating similar thickness-dependent effect.

Density-based topology optimization has emerged as a robust and versatile technique for designing structures with optimal material distribution. Conventional solid isotropic material with penalization (SIMP) approaches (Bendsøe and Sigmund 2004) penalize intermediate densities to promote a near-binary material distribution. Furthermore, a Heaviside projection filter (Guest et al. 2004) is incorporated to introduce a characteristic length scale in the design space, enabling feature size control and mitigating mesh dependency issues.

Recently, coating filters have been introduced in the literature to model surface layer effects in optimized structures. Clausen et al. (2015) presented a topology optimization framework for structures with stiff coating. Subsequent to this investigation, multiple adaptations of coating filters have been explored (Groen et al. 2019; Luo et al. 2019; Yoon and Yi 2019; Suresh et al. 2020). Wang and Kang (2018) further extended these concepts by introducing a level-set method to model coated structures within a shape and topology optimization framework. While these methods can incorporate thickness-dependent effects, they were not specifically developed for this purpose. In addition, Tuna and Trovalusci (2022) conducted topology optimization of two-dimensional plates, incorporating Eringen’s nonlocal theory of elasticity (Eringen 1987), to capture size effects. However, unlike traditional nonlocal elasticity models that compute a nonlocal stress measure at each point in the spatial domain, our approach employs a simpler nonlocal formulation for the material’s elastic constant, enabling the computation of a spatially varying elastic modulus at the surface.

This study aims to develop a density-based topology optimization framework tailored for thin-walled structures, incorporating the thickness-dependent nature of material properties. The nonlocal integral model from Jaksch et al. (2022) is employed, which effectively captures this effect by inducing surface grading, motivated by the non-uniform heat distribution in the laser-sintering process. This enables the direct integration of material property variations into topology optimization, enhancing manufacturability while preserving structural efficiency. The framework is particularly well suited for optimizing lattice microstructures (Wu et al. 2021), where periodic cell dimensions typically range from 5 to 20 mm.

The remainder of this paper is organized as follows. Section 2 presents the problem formulation. This includes the definition of the nonlocal material model and the optimization problem. Section 3 presents the results of the numerical experiments. Finally, Section 4 concludes the paper.

2 Problem formulation

This section presents the nonlocal material model and optimization framework used in this study. The nonlocal integral model for Young’s modulus variation follows the approach of Jaksch et al. (2022), which captures thickness-dependent stiffness degradation through a spatially weighted integral. However, that work focused solely on material modeling without integrating it into topology optimization. Our contribution is to extend this model into a density-based topology optimization framework, allowing it to directly influence material distribution.

2.1 Material model

Consider $Ω \in R^{d}, d \in {2, 3}$ as a homogeneous, bounded physical domain. Let $I_{Ω}$ be an indicator function defining the interior of the domain, given by

(1)

Given a reference Young’s modulus $E_{0}$ and a kernel function $ψ : R^{d} \to R$ , the Young’s modulus at a point $x x \in Ω$ is defined as

(2)

where $ϕ (x x)$ represents the nonlocal term, expressed as the convolution integral:

(3)

This nonlocal term introduces a dependency on the spatial distribution of Young’s modulus. The parameter $ζ \in [0, 1]$ serves as a fraction coefficient, offering control over the surface grading profile. Particularly, for $ζ = 1$ , the material exhibits a uniform Young’s modulus distribution. The kernel function $ψ$ captures the nonlocal effect, showcasing a decaying nature that vanishes beyond a certain limit. It satisfies the normalization condition:

(4)

In this study, we adopt the following expression for the kernel function:

(5)

which is widely used in density-based topology optimization literature. Here, $δ > 0$ defines the region of influence of the kernel, determining the extent to which nonlocal interactions impact the material properties.

To elucidate the impact of the selected material model on the effective Young’s modulus concerning wall thickness, we analyze a one-dimensional specimen with thickness denoted by t, as illustrated in Fig. 1. The homogenized effective Young’s modulus at a point $x \in [- t / 2, t / 2]$ , denoted as $E_{eff}$ , is computed as the average of the varying Young’s modulus E across the thickness of the specimen, expressed as

(6)

Figure 2 shows the dependency of the effective Young’s modulus on the thickness of the specimen, considering $δ = 0.4 mm$ (Jaksch et al. 2022) and varying values of the fraction coefficient $ζ$ . The integrals in Eqs. (3), (5), and (6) are computed numerically using the mid-point method. For values of $ζ < 1$ , the effective Young’s modulus significantly decreases for thickness below 1 mm. For $ζ = 1$ (fully local behavior), the effective Young’s modulus is independent of the specimen’s thickness.

Figure 3 depicts the spatial variation of Young’s modulus across the thickness of the one-dimensional specimen, highlighting distinct profiles associated with varying thickness values. For specimens with a thickness $t \geq 2 δ$ , the graded surface layer maintains a constant thickness of $δ$ . In this scenario, the Young’s modulus at the surface is $0.5 E_{0}$ and increases gradually until reaching its maximum value of $E_{0}$ within the graded layer. However, for specimen thicknesses below $2 δ$ , the graded surface layer encompasses the entire thickness, leading to a decrease in both the maximum and minimum values of Young’s modulus with diminishing thickness.

2.2 Topology optimization

The goal of this study is to design structures with graded surfaces in terms of elastic modulus, while ensuring a near-binary material distribution. The following section outlines the optimization framework.

2.2.1 Design parameterization

The design domain is discretized using N finite elements, followed by the application of a density-based topology optimization method within the spatially discrete framework. In conventional topology optimization using the SIMP method, a piecewise-constant, scalar density field, represented by $ρ_{e}, e = 1, \dots, N$ , is employed. Here, $ρ_{e} = 0$ denotes a void element, $ρ_{e} = 1$ indicates a solid element, and intermediate values of $ρ_{e} \in (0, 1)$ represent gray regions. These intermediate regions are penalized in the optimization process, promoting a near-binary material distribution. These values are termed “densities” as they are combined with elemental volumes $v_{e}, e = 1, \dots, N$ , to give the total mass of the structure.

However, a drawback of this penalization approach is the emergence of checkerboard pattern in the optimized design when piecewise-linear finite element interpolation is used for the displacement field. To overcome this issue, a pseudo-density field, given by $μ_{e}, e = 1, \dots, N$ , is introduced. The pseudo-densities act as control variables in the optimization process, and lack a direct physical interpretation. A convolution filter is applied to the pseudo-density field which effectively mitigates the occurrence of checkerboard pattern (Bruns and Tortorelli 1998).

The introduction of the filter also establishes a length scale in the optimized design, addressing the issue of the optimizer converging to designs dependent on the resolution of the finite element discretization. However, if the filter’s radius is excessively large, it results in significant gray regions. To overcome this, the filtered values, denoted by ${\tilde{μ}}_{e}, e = 1, \dots, N$ , are translated to a near-binary density field using a smooth approximation of the Heaviside step function.

In this study, we employ this strategy, utilizing the commonly employed density filter on pseudo-densities $μ_{e}, e = 1, \dots, N$ , viz.,

(7)

where filter weights $w_{e i}$ are computed based on the distances between the finite element centroids $x x_{e}, e = 1, \dots, N$ , as

(8)

with R representing the filter radius. Subsequently, the intermediate values ${\tilde{μ}}_{e}, e = 1, \dots, N$ , undergo projection using a smooth approximation of the Heaviside step function to give material densities $ρ_{e}, e = 1, \dots, N$ , expressed as

(9)

The parameter $β > 0$ governs the sharpness of the projection, while $η \in [0, 1]$ acts as the threshold parameter. Notably, when $η$ assumes values of 0 and 1, the projection corresponds to the Heaviside step filter (Guest et al. 2004) and modified Heaviside filter (Sigmund 2007), respectively. In the limit where $β \to \infty$ , these threshold values enable control over the length scale for the solid and the void phases, respectively.

2.2.2 Surface grading of Young’s modulus

To incorporate surface grading of Young’s modulus, we apply the nonlocal material model outlined in Sect. 2.1 to the density field $ρ$ . For the approximation of the integral in Eq. (3), we utilize mid-point integration, specifically one-point Gaussian quadrature. Assuming the support of the kernel function $ψ (\cdot; x x_{e})$ can be contained within the design domain for atleast one $x x_{e}, e = 1, \dots, N$ , the resulting integrated density ${\bar{ρ}}_{e} \in [0, 1]$ approximating $ϕ$ is given by

(10)

Here, the $max$ operator in the denominator ensures that the nonlocal model behaves in a physically consistent manner at the edges of the domain.

Since ${\bar{ρ}}_{e}$ can be greater than zero even for void elements ( $ρ_{e} = 0$ ), we introduce a scalar field represented by $α_{e} \in [0, 1], e = 1, \dots, N$ , where

(11)

The parameter $α_{e}$ acts as a scaling coefficient, which adjusts the reference Young’s modulus $E_{0}$ to obtain the elemental Young’s modulus. The inclusion of the exponent $p > 1$ aims to penalize intermediate densities $ρ_{e}$ , and encourages the emergence of a near-binary material distribution. In conventional compliance minimization problems, a common choice is $p = 3$ as used in the literature (Bendsøe and Sigmund 2004). However, due to the multiplication of the term $ρ_{e}^{p}$ by a linear function of the filtered density ${\bar{ρ}}_{e}$ in the above expression, lower values of p depending on $ζ$ are more appropriate to prevent excessive penalization that potentially slows down convergence in the optimization procedure.

The elemental Young’s modulus corresponding to finite element e is then expressed as

(12)

where a minimal stiffness coefficient of $κ = 10^{- 9}$ is incorporated for void elements ( $ρ_{e} = α_{e} = 0$ ) to circumvent singularity issues in the global stiffness matrix during the numerical implementation.

2.2.3 Optimization problem

To exemplify the material model within the context of a structural design problem, we examine the standard compliance minimization problem with a structural volume constraint. The optimization problem is formulated as follows:

(13)

Here, $c : R^{N} \to R$ represents the global compliance, encompassing the overall stiffness of the structure. The vectors $u \in R^{M}$ and $f \in R^{M}$ denote the nodal displacements and forces, respectively. The structure’s volume is defined by the function $g : R^{N} \to R_{\geq 0}$ , and $V > 0$ is the prescribed volume. For a given design $μ μ$ and external load $f$ , $u = u (α α (μ μ))$ is the solution of the state problem:

(14)

Here, $K_{e}$ denotes the elemental stiffness matrix associated with element e, assuming a unit Young’s modulus.

For solving the topology optimization problem (Eq. 13), a solver based on the method of moving asymptotes (MMA) (Svanberg 1987) is used, implemented in the ParOpt (Chin et al. 2019) library. The gradient of the objective function as required by the optimizer is computed using the adjoint method (Bendsøe and Sigmund 2004).

The sensitivity of the compliance c with respect to the design parameter $μ_{e}$ is expressed as

(15)

where the derivative $\partial c / \partial α_{e}$ using the adjoint method is formulated as

(16)

The other partial derivatives for $i, j = 1, \dots, N$ are given as

(17)

(18)

(19)

(20)

where $δ_{i j}$ is Kronecker delta. Similarly, using the chain rule, the derivative $\partial g / \partial μ_{e}$ is given by

(21)

3 Numerical examples

To evaluate the effectiveness of the proposed topology optimization framework, we apply it to three case studies: a 2D cantilever beam, a 2D bridge structure, and a 3D pressure box. These examples were chosen because, without length scale control, they exhibit thin features after topology optimization. Our nonlocal integral model naturally suppresses such features, ensuring manufacturable designs while maintaining structural performance. A plane-strain condition is assumed. We set the reference Young’s modulus to $E_{0} = 1$ MPa, and the Poisson’s ratio is assumed to be 0.3.

The geometrical dimensions of the two structures are illustrated in Fig. 4. The domain is discretized using bilinear quadrilateral elements with thickness of 1 mm. The cantilever beam is fixed at the left end, while a static traction of 1 MPa distributed over 1 mm length is applied at the right end. For the bridge structure, a distributed load of 1 MPa is applied on the top surface, while the bottom corners are held fixed. The regions in dark gray indicate non-design regions, introduced to have full elastic stiffness at the support and loading points. The volume fraction for the design domain is set to 0.4 for the cantilever beam problem and 0.2 for the bridge problem.

In all the test cases, we use $ζ = 0$ , that corresponds to fully nonlocal behavior. The case with $ζ = 1$ corresponds to pure local behavior which is the same as using $δ = 0$ .

3.1 Size effect

First, we examine the influence of the overall design domain size on the optimized designs and their corresponding structural compliance. Optimizations are conducted on distinct domain sizes, scaled by factors of $s = 1, 2, 3$ and 4. These correspond to domain dimensions of $15 \times 10$ ${mm}^{2}, 30 \times 20$ ${mm}^{2}, 45 \times 30$ ${mm}^{2},$ and $60 \times 40$ ${mm}^{2}$ for the cantilever beam. For the bridge structure, the domain sizes are $20 \times 10$ ${mm}^{2}, 40 \times 20$ ${mm}^{2}, 60 \times 30$ ${mm}^{2},$ and $80 \times 40$ ${mm}^{2}$ , respectively.

Despite variations in domain size, the finite element discretization is consistently executed with the same number of elements. Consequently, the mesh size h is proportionally scaled by the factor s. For the coarsest mesh in both the cantilever beam and bridge cases, $h = H = 1.5625 e - 2 mm$ . The radius of the Heaviside projection filter is set to $R = 10 h$ . The continuation parameter $β$ is systematically increased at intervals of every 50 optimization iterations, taking values of $β = 1, 2, 4, 8, 16, 32$ , and 64. After reaching the highest value of $β$ , the optimization process is continued until convergence, with a maximum of 1000 iterations. The threshold parameter $η$ is chosen to be 0.5.

For optimization, two distinct values of the kernel parameter $δ$ are considered. The first value, $δ = 0$ , corresponds to the standard SIMP approach with the Heaviside projection filter. Throughout the optimization, the penalization constant p remains fixed at $p = 3$ . The second value, $δ = 0.4 mm$ , enables the emulation of surface grading effects. For this case, the penalization constant is selected as $p = 2.5$ .

Figure 5 showcases the optimized designs. The Young’s modulus distribution and compliance values depicted are derived from a postprocessing stage involving volume-preserving thresholding followed by the employment of the nonlocal material model with $δ = 0.4 mm$ . It is noteworthy that the designs incorporating surface grading exhibit slightly lower compliance compared to those generated without the explicit consideration of surface grading in the optimization procedure. Moreover, as the domain size increases, the influence of surface grading diminishes, leading to nearly identical designs for both $δ$ values. Figure 6 shows the relationship between domain scale and compliance, with compliance values normalized by those obtained for designs optimized and postprocessed with $δ = 0$ . This normalization enables direct comparison across different domain sizes. The results show that normalized compliance increases as domain size decreases, indicating greater sensitivity to surface grading in smaller-scale structures. Additionally, while the compliance difference between optimizations with ( $δ = 0.4 mm$ ) and without ( $δ = 0$ ) surface grading remains small, it grows as domain size decreases, emphasizing the method’s role in reducing excessively thin features and enhancing manufacturability.

The horizontal dashed line in Fig. 6 represents the compliance trend for the standard SIMP approach. Since we use a mesh-dependent density filter, the optimized designs remain independent of domain size. This means that, under standard SIMP, scaling the domain does not influence compliance trends, unlike in the proposed nonlocal approach, where surface effects become more pronounced in smaller structures.

This trend is expected. As domain size increases while keeping the kernel radius $δ$ and discretization fixed, the relative effect of surface grading decreases. Larger structures naturally exhibit thicker features, reducing the proportion affected by stiffness degradation near surfaces. In contrast, smaller domains contain more thin features, making surface grading effects more pronounced. This aligns with theoretical expectations, as surface effects are generally more influential in fine-scale structures than in large-scale topologies.

3.2 Mesh dependency

Next, we study optimized designs without incorporating the Heaviside projection filter, and the necessity for a continuation scheme for $β$ in the optimization procedure. While this method yields designs that exhibit sensitivity to mesh variations and lack integration of feature size control, it concurrently offers the advantage of a simplified algorithm. By setting R to a small yet sufficiently large value (e.g., $R = 1.3 h$ ), checkerboarding issues are effectively mitigated

Figures 7 and 8 visually present the optimized designs for the cantilever beam and bridge, respectively. These designs are generated under varying degrees of mesh refinement and are juxtaposed against the standard SIMP method, corresponding to $δ = 0$ . Numerical experiments affirm that, with moderately refined meshes, the resultant designs exhibit reduced variability when compared to those generated by the standard SIMP method.

However, there are two main downsides when using highly refined meshes. First, there is a higher chance of getting stuck in local minima, making optimization more challenging. Second, as the mesh gets finer, we start to see more complex structures (rank-2 microstructures) in the designs. This complexity adds challenges to the optimization process (Sigmund and Petersson 1998).

3.3 Extension to 3D

To demonstrate the applicability of the proposed framework in three dimensions, we examine its application to the optimization of reinforcement structure for a symmetrical box subjected to an internal pressure. The inner face has an edge length of 10 mm and shell thickness of $0.5 mm$ , as shown in Fig. 9. Exploiting the geometrical symmetry, computational efficiency is gained by simulating and optimizing only one-eighth of the box. A kernel radius of $δ = 0.4 mm$ is selected, and a mesh size of $h = 0.0625 mm$ is employed, yielding a ratio of $δ / h = 6.4$ . The optimization problem is solved using 2473984 design variables. The filter radius is chosen as $R = 1.3 h$ and no Heaviside projection is employed during the optimization process. We assume a fraction coefficient of $ζ = 0$ and constrain the volume of the full structure to 800 ${mm}^{3}$ .

Figures 10 and 11 show the optimized design, which appears in the form of a spherical structure with thin members connecting the spherical shell and the box.

4 Conclusion

In this work, we have presented a framework for the topology optimization of structures featuring graded surfaces, leveraging a nonlocal methodology that captures thickness-dependent stiffness variations. This approach addresses a critical aspect of additively manufactured structures by incorporating material property grading directly into the optimization process.

We have demonstrated how this framework influences the distribution of Young’s modulus across the thickness in one- and two-dimensional structures, revealing nuanced variations in material properties. The key advantage of our approach is its ability to optimize structures while preserving graded surface effects, which is not readily achievable with conventional SIMP-based topology optimization methods.

Our findings suggest that incorporating material property variations into the optimization process leads to more structurally realistic designs, particularly in small-scale, thin-walled structures, where thickness-dependent effects are significant. While compliance improvements may be moderate, the trend toward enhanced stiffness and reduced deformation aligns with theoretical expectations, demonstrating the benefit of integrating nonlocal modeling in topology optimization. Additionally, our method also works without explicit feature size control, as surface degradation effects naturally suppress the formation of thin features. However, if a specific minimum feature size is required, additional length scale control can still be incorporated.

A limitation of the proposed framework is the computational cost of nonlocal integration, particularly in distributed computing. The current implementation incurs significant overhead due to neighbor search and storage, with run-time and memory allocation scaling with the number of neighbors within the integration region. A more scalable alternative is the use of Helmholtz-type PDE-based filter (Lazarov and Sigmund 2011), which provides a computationally efficient approach for enforcing nonlocal effects. Developing and integrating a PDE-based nonlocal model into topology optimization remains an important direction for future research.

This framework can be extended in several directions. Future work could integrate material anisotropy, as explored in Suresh et al. (2020), that could further enhance the versatility and applicability of the approach. Another promising direction is the extension to failure-driven topology optimization. A nonlocal damage model (Peerlings et al. 1996) could be incorporated to simulate strength degradation in thin-walled structures, where a separate nonlocal kernel function is calibrated using experimental strength data. Damage constraints could then be integrated into the topology optimization framework, as explored in Russ and Paulino (2023), providing a more comprehensive optimization approach that accounts for both stiffness and failure resistance in additively manufactured structures.

Related Articles

Contact us