Topology Optimization for Digital Light Projector Additive Manufacturing Addressing the In-Situ Structural Strength Issue

Abstract

A topology optimization approach is proposed for the design of self-supporting structures for digital light projector (DLP) 3D printing. This method accounts for the adhesion forces between the print part and the resin base during DLP printing to avoid failure of the part due to stress concentration and weak connections. Specifically, the effect of the process-related adhesion forces is first simulated by developing a design variable-interpolated finite element model to capture the intricate mechanical behavior during DLP 3D printing. Guided by the process model, a stress-constrained topology optimization algorithm is formulated with both the SIMP and RAMP interpolation schemes. The interpolations on the stress term and the design-dependent adhesion load are carefully investigated. A sensitivity result on the P-norm stress constraint is fully developed. Finally, the approach is applied to several 2D benchmark examples to validate its efficacy in controlling the process-caused peak P-norm stresses. The effects of alternating between the SIMP and RAMP interpolations and changing the stress upper limits are carefully explored during the numerical trials. Moreover, 3D printing tests are performed to validate the improvement in printability when involving the process-related P-norm stress constraint.

1. Introduction

Additive manufacturing reverses the conventional subtractive manufacturing process that slices the target part into layers and performs the manufacturing in a bottom-up layer-by-layer approach. The benefit of this novel approach is multifaceted: (i) complex-shaped parts can be manufactured similarly to their simple-shaped counterparts that decouple the proportional relationship between the manufacturing difficulty and the geometric complexity [1]; (ii) parts with tremendous interior cavities such as lattice structures that cannot be handled by conventional manufacturing are enabled for physical fabrication with additive manufacturing [2]; (iii) tedious tooling development (e.g., molds and fixtures) can be eliminated, which reduces the manufacturing preparation, saving both preparation time and costs; (iv) a variety of additive manufacturing processes have been developed and matured for industrial production that cover nearly all types of materials, including polymers, ceramics, metals and composites. The ISO/ASTM 52900:2015 has standardized the terminologies and principles of the diversified processes [3].

On the other hand, additive manufacturing is still a relatively expensive process compared to the conventional approaches due to the higher demands for raw materials and longer production cycles [4]. Hence, the success of applying additive manufacturing often leads to producing novel products dedicated to the layered manufacturing approach, i.e., equipped with extraordinary functionality or light weight by exploring the greatly expanded geometric design space that qualifies the costly additive production [4,5]. For the same reason, topology optimization is tightly connected to additive manufacturing, relying on algorithms and computational intelligence to achieve designs for additive manufacturing. The design results from topology optimization in general demonstrate excellent physical performance while utilizing only a limited quantity of materials, ensuring a light weight. Therefore, topology optimization for additive manufacturing has been of great significance in the past decade and interested readers can refer to a few dedicated review papers [6,7] for details.

Topology optimization for additive manufacturing takes into account the process characteristics to enhance the printability of the topological designs. The printability here can indicate a higher manufacturing success rate, lower manufacturing time or costs, a reduced gap between the as-designed and as-fabricated structural performance, etc. For instance, topology optimization for metal additive manufacturing takes the residual stress and distortion constraints into account [8,9,10,11], relieving the residual stress concentration to prevent cracks and reducing the z-axis warpage to avoid collisions to the powder brush (for powder bed fusion processes). Topology optimization for filament-based 3D printing processes, especially fused deposition modeling (FDM) for fiber-reinforced composites, includes both the topological variables and the path direction variables as the optimization targets [12,13,14], performing the concurrent optimization of the design with material anisotropies [15], which closes the gap between design and manufacture. Overhang-free topology optimization [16,17,18,19,20] controls the inclination angles of the overhanging features to eliminate the need for support structures, reducing the support producing time and simultaneously the cost [21,22]. The overhang-free optimization is generally applicable to all types of support-involved additive manufacturing processes. More issues, like the multi-scale design problem [23,24,25,26,27] and the post-machinability problem [28,29,30], have also been addressed.

Digital light projector (DLP), polymerizing the photosensitivity resins through the layered projection of UV light following the sliced 2D images, is a type of vat polymerization 3D printing process. The application scope of the DLP process is quite widespread since it is able to cure polymers, ceramics, and even metals with an excellent curing resolution [3]. The application cases for aerospace, aeronautics, and biomedical structures are widely reported. However, topology optimization for DLP 3D printing seems to be less discussed. Practically, the topology optimization results can be directly printed with the DLP machine and only sparse supports are required for process safety compared to other 3D printing techniques. On the other hand, process failure in DLP printing is frequently encountered with the repeatedly applied adhesion forces. Referring to Figure 1, the polymerized resin is bonded to the build bed but also glued to the tank base and, hence, adhesion forces take effect on the polymerized resin when separating it from the tank base to prepare the next layer for UV curing. Hence, the adhesion forces are applied to all layers and are likely to damage the part if weak connections exist. In summary, topology optimization for DLP 3D printing addressing the in-process structural strength issue deserves careful investigation.

To ensure the structural strength, stress-constrained topology optimization is generally performed. The anisotropic failure criteria [31,32] are applied to strength-based topology optimization due to the process-induced material anisotropy. Given the algorithm details under the density-based framework, the qp method [33,34] is widely adopted to prevent the stress singularity problem, i.e., eliminating high stress levels in the low-density elements. The P-norm or Kreisselmeier–Steinhauser function is applied to aggregate the many local stress constraints into one global peak stress constraint [35,36,37,38], and iterative adjustment to the P-norm approximation is also necessary to close the gap between the approximated and the exact peak stresses [39]. The other numerical issue is the high nonlinearity of the optimization problem, where the stress levels can be sharply perturbed by minor material distribution changes [35]. Hence, all parameters involved in configuring and solving the optimization problem should be fine-tuned. If performed under other topology optimization frameworks, the algorithm details would differ significantly. For instance, the black and white structural representation with the level set method avoids addressing the stress singularity issue [40,41,42]. The sharp solid-void transition for the Bidirectional Evolutionary Structural Optimization (BESO) method also eliminates this issue [43].

In this research, the stress-constrained topology optimization is configured differently from the majority of existing works. Specifically, apart from the work load, the process-induced load is applied in a layer-by-layer manner, indicating that, for every intermediate state, there is a coupled stress field to restrict the peak magnitude. Hence, multiple P-norm approximations on the peak stress are applied, each corresponding to a 3D printing intermediate state. Additionally, the process-induced load is design-dependent since more materials to be cured in the current layer will lead to a larger total force dragging the intermediate part. The stress-related problems involving design-dependent loads (pressure load or self-weight) are even more sophisticated [44,45,46], and Rational Approximation of Material Properties (RAMP) interpolation [47] is often used for this type of problem [8]. Hence, both the Solid Isotropic Material with Penalization (SIMP) and RAMP interpolations are investigated. More algorithm details are given in the following sections.

2. Topology Optimization Considerations for DLP 3D Printing

In the DLP 3D printing process, as shown in Figure 1, due to the upward motion of the print bed, adhesive forces are generated when separating the in-process workpiece from the tank base. Fractures or even complete failure may occur if weak structural connections exist, which should be avoided. Material design before printing can alleviate this issue [48], and, alternatively, topology optimization involving process-related stress constraints may completely eliminate this problem.

2.1. DLP Process Model

As illustrated in Figure 2, the domain $\mathsf{\Omega }$ is divided into m layers, with each layer defined as ${\Omega }_i$, where 1 ≤ $i$ ≤ m. The bottom is fixed on the print bed, and the material curing process proceeds layer-by-layer in the bottom-up direction. Hence, the domain $\Omega$ is divided into three subdomains, ${\Omega }_d$, ${\Omega }_e$, and ${\Omega }_h$, where ${\Omega }_d$ is the collection of layers that have already been printed at the current stage, ${\Omega }_e$ is the print layer that bears the adhesion forces in the current load step, and ${\Omega }_h$ is the collection of layers that have not yet been printed.

During the printing process, the force-bearing layer ${\Omega }_e$ is subjected to adhesion forces due to the lifting of the print platform and the detachment from the vat base. The adhesion forces act only on the polymerized materials and, hence, the force distribution is related to the design variables, serving as a design-dependent load and interpolated with the SIMP model specifically as

$$\left\{\begin{array}{cc}{\delta }_e^i=1\hfill & when {\Omega }_i\in {\Omega }_e\\ {\delta }_e^i=0\hfill & otherwise\end{array}\right.$$

(1)

$$f_{element}={\delta }_e^i{\rho }_e^{3}ga$$

(2)

where ${\delta }_e^i$ is the indicator function representing domain ${\Omega }_e$, $f_{element}$ is the total adhesion force applied to an element, $g$ is the adhesion coefficient [49,50], and $a$ is the length of the rectangular cell.

In topology optimization, equivalent nodal forces are required to represent the adhesion forces with the finite element method. Here, four-node quadrilateral elements are used for the finite element analysis and the equivalent nodal force at the j-th node of a four-node quadrilateral element can be obtained as

$$f_{iy}^e=\int {\delta }_e^i{N}_j{\rho }_e^{3}gadA={\delta }_e^i{\rho }_e^{3}ga\int N_{j}dA=\frac{1}{4}{\delta }_e^i{\rho }_e^{3}g{a}^{2}b$$

(3)

where $N_j$ is the shape function of the j-th node of the four-node rectangular element;${\rho }_e$ is the physical density of the element; $b$ is the width of the rectangular cell.

The progressive material curing process is simulated by configuring the element activation strategy in the finite element model. Specifically, ${\Omega }_h$ is filled with substitute materials and will be activated sequentially in the bottom-up manner with the following indicator variable, ${\theta }_e^i$, which is the indicator function representing domain ${\Omega }_d$.

$$\left\{\begin{array}{cc}{\theta }_e^i=1\hfill & when {\Omega }_i\in {\Omega }_d\\ {\theta }_e^i=0\hfill & otherwise\end{array}\right.$$

(4)

When printing the i-th layer, the global stiffness matrix is the sum of the element stiffness matrices in ${\Omega }_d$, which can be explicitly expressed using the indicator variable:

$${\mathit{K}}_{\mathit{i}}={\displaystyle \sum }_{e=1}^{nele}({\theta }_e^i{\mathit{L}}_{\mathit{e}}^{\mathit{T}}{\mathit{K}}_{\mathit{e}}{\mathit{L}}_{\mathit{e}})$$

(5)

In this context, the matrix ${\mathit{L}}_{\mathit{e}}$ collects the nodal displacements $u_e$ of the e-th element from the global displacement vector $\mathit{U}$, satisfying the equation ${\mathit{u}}_{\mathit{e}}={\mathit{L}}_{\mathit{e}}\mathit{U}$. ${\mathit{K}}_{\mathit{e}}$ represents the stiffness matrix of the element, and its calculation is conducted as follows:

$${\mathit{K}}_{\mathit{e}}=\int {\mathit{B}}^{\mathit{T}}{\mathit{D}}_{\mathit{e}}\mathit{B}d{\Omega }_e={\eta }_e^E\int {\mathit{B}}^{\mathit{T}}{\mathit{D}}_0\mathit{B}d{\Omega }_e$$

(6)

Here, ${\eta }_e^E$ represents the interpolation function for the elastic modulus, $\mathit{B}$ denotes the strain–displacement matrix, and ${\mathit{D}}_0$ corresponds to the constitutive matrix of the solid material. The displacement ${\mathit{U}}_{\mathit{i}}$ due to the adhesion load at printing stage i is calculated as follows:

$${\mathit{K}}_{\mathit{i}}{\mathit{U}}_{\mathit{i}}={\mathit{F}}_{\mathit{i}}$$

(7)

In the implementation, multiple sliced layers are consolidated to eliminate the need for an extremely fine mesh and reduce the number of calls on finite element analysis.

2.2. AM Filter for Fabricable Self-Support Structures

A schematic diagram of the amplitude filter is provided in Figure 3 to illustrate its working principle. In this filter, the cells located in the lower layer with indices (i − 1, j − 1), (i, j − 1), and (i + 1, j − 1) are defined as the support region for cell (i, j). If there is no material in the support region, the amplitude filter will remove material from cell (i, j); otherwise, material is allowed to remain in cell (i, j). Through this rule, the amplitude filter achieves a 45° overhang self-supporting characteristic. The mathematical representation of the amplitude filter is shown below:

$${\xi }_{i,j}=\left\{\begin{array}{c}\mathit{min}\left({\overline{\rho }}_{i,j}, \mathit{max}\left({\overline{\rho }}_{i-1,j-1},{\overline{\rho }}_{i,j-1},{\overline{\rho }}_{i+1,j-1}\right)\right) \left(i,j\right)\in {\Omega }_u\\ {\overline{\rho }}_{i,j} \left(i,j\right)\in {\Omega }_b\end{array}\right.$$

(8)

The projection field $\overline{\rho }$ has its differential form expressed as follows:

$${\xi }_{i,j}=smin\left({\overline{\rho }}_{i,j}, smax\left({\overline{\rho }}_{i-1,j-1},{\overline{\rho }}_{i,j-1},{\overline{\rho }}_{i+1,j-1}\right)\right) \left(i,j\right)\in {\Omega }_u$$

(9)

The operator “$smin$“ is defined as follows:

$$smin\left(a,b\right)=\frac{1}{2}\left(a+b-\sqrt{\left({\left(a-b\right)}^2+{ϵ}_s\right)}+{ϵ}_s{}^2\right) \left(i,j\right)\in {\Omega }_u$$

(10)

The operator “$smax$” represents the maximum function “P-Qmax”, which is used to calculate the maximum value of elements within the support region.

$$smax\left(a,b,c\right)=\sqrt[Q]{a^P+b^P+c^P} \left(i,j\right)\in {\Omega }_u$$

(11)

The parameter ${ϵ}_s=1\times {10}^{-4}$ controls the approximation accuracy. The usage of $P=40$ and $Q=P+\frac{\log 3}{\log \frac{1}{2}}$ is followed.

3. Interpolation Scheme and Total Stress Calculation

Density-based topology optimization conventionally uses interpolation schemes like Solid Isotropic Material with Penalization (SIMP) and Rational Approximation of Material Properties (RAMP) to relate the material properties to the element densities (topological variable). This study employs both the SIMP and RAMP interpolations to gain insights into their differences and similarities for adhesion-aware topology optimization. SIMP uses a power-law function to penalize intermediate densities, while RAMP uses a rational function to provide a smooth transition between solid and void phases. We implement both schemes in this work to compare their performance in mitigating adhesion-induced stress concentration, a key challenge in DLP-type additive manufacturing. The results will disclose the benefits and limitations of both interpolations in handling in-situ structural strength constraints in topology optimization for polymer structures made by DLP 3D printing. The specific interpolation schemes are given below:

$${\mathrm{E}}_e={\mathrm{E}}_{min}+{\eta }_e^E\left({\mathrm{E}}_0-{\mathrm{E}}_{min}\right)$$

(12)

SIMP: ${\eta }_e^E={\rho }_e{}^{p\_SIMP}$

RAMP: ${\eta }_e^E=\frac{{\rho }_e}{1+\left(1-{\rho }_e\right)p\_RAMP}$

${\mathrm{E}}_e$ is the interpolated elastic modulus for element e; $p\_SIMP$ is the penalty parameter of the SIMP interpolation; $p\_RAMP$ is the RAMP interpolation penalty parameter; ${\rho }_e$ is the physical density of element e; ${\mathrm{E}}_0$ represents the elastic modulus of the solid material; ${\mathrm{E}}_{min}$ is the elastic modulus for voids.

The stress vector ${\mathit{\sigma }}_{\mathit{e}}$ for element e is calculated through

$${\mathit{\sigma }}_{\mathit{e}}={\theta }_e^i{\mathit{D}}_0{\mathit{B}}_{\mathit{e}}{\mathit{u}}_{\mathit{e}}$$

(13)

where ${\mathit{D}}_0$ is the stiffness matrix of the solid material; ${\mathit{B}}_{\mathit{e}}$ is the strain–displacement matrix; ${\mathit{u}}_{\mathit{e}}$ is the local displacement vector. The two-dimensional stress vector ${\mathit{\sigma }}_{\mathit{e}}$ is

$${\mathit{\sigma }}_{\mathit{e}}={\left({\sigma }_{ex},{\sigma }_{ey},{\tau }_{exy}\right)}^T$$

(14)

The relaxed stress measure ${\widehat{\mathit{\sigma }}}_{\mathit{e}}\left({\rho }_e\right)$ is expressed as

$${\widehat{\mathit{\sigma }}}_{\mathit{e}}\left({\rho }_e\right)=\eta \left({\rho }_e\right){\mathit{\sigma }}_{\mathit{e}}$$

(15)

A general stress punishment scheme is given as

$$\eta \left({\rho }_e\right)={\left({\rho }_e\right)}^q$$

(16)

where $q$ is a non-negative stress relaxation parameter. For full material cells, ${\widehat{\mathit{\sigma }}}_{\mathit{e}}$ is equal to ${\mathit{\sigma }}_{\mathit{e}}$, and for an empty cell, it satisfies

$$\underset{{\rho }_e\to 0}{\lim }{\widehat{\mathbf{\sigma }}}_{\mathbf{e}}\left({\rho }_e\right)=0$$

(17)

The two-dimensional von Mises stress is defined by

$${\widehat{\sigma }}_{vm,e}={\left({\sigma }_{ex}{}^2+{\sigma }_{ey}{}^2-{\sigma }_{ex}{\sigma }_{ey}+3{\tau }_{exy}^2\right)}^{\frac{1}{2}}$$

(18)

A standard P-norm global stress measurement is used to approximate the maximum stress as follows:

$${\sigma }_{PN}={\left({{\displaystyle \sum }}_{e=1}^{nele}{\widehat{\sigma }}_{vm,e}{}^p\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$p$}\right.}$$

(19)

where ${\widehat{\sigma }}_{vm,e}$ is the von Mises stress at the centroid of element e, and $p$ is the P-norm aggregation parameter. As $p\to \infty$, the P-norm approaches the maximum value of ${\widehat{\sigma }}_{vm}$.

$$\mathrm{Max}\left({\widehat{\sigma }}_{vm}\right)\le {\left({{\displaystyle \sum }}_{e=1}^{nele}{\widehat{\sigma }}_{vm,e}{}^p\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$p$}\right.}$$

(20)

In general, a larger $p$ provides a more accurate approximation of the maximum von Mises stress. However, if $p$ is excessively large, it can lead to an ill-posed problem and cause severe oscillations during the optimization process. Therefore, it is important to select an appropriate $p$ value to ensure smooth convergence while sufficiently approximating the maximum stress.

4. Formulation and Solution of Optimization Problems

4.1. Optimization Problem Formulation

The objective function of the optimization problem is to minimize the compliance of the structure under the maximum stress constraint and volume constraint, and the mathematical expression of the problem is as follows:

$$\left\{\begin{array}{c}find: \mathit{\rho }\\ minimize:C={\mathit{U}}_{\mathit{c}}^{\mathit{T}}{\mathit{K}}_{\mathit{c}}{\mathit{U}}_{\mathit{c}}\\ subject to:\left\{\begin{array}{c}{\mathit{K}}_{\mathit{i}}{\mathit{U}}_{\mathit{i}}={\mathit{F}}_{\mathit{i}}\left(i=1,2,\dots ,m\right)\\ {\mathit{K}}_{\mathit{c}}{\mathit{U}}_{\mathit{c}}={\mathit{F}}_{\mathit{c}}\\ \frac{V\left(\mathit{\rho }\right)}{V_0}=f\\ {\sigma }_{PN}^i\le \sigma \\ 0<{\rho }_{min}\le \forall {\rho }_e\le 1\end{array}\right.\end{array}\right.$$

In this equation, ${\mathit{U}}_{\mathit{c}}$ represents the nodal displacement field induced by mechanical external loads, ${\mathit{K}}_{\mathit{c}}$ denotes the overall stiffness matrix of the structure, and $V\left(\rho \right)={{\displaystyle \sum }}_{e=1}^{nele}{\stackrel{=}{\rho }}_e$ represents the material volume. $V_0$ denotes the total volume of the design domain. Setting ${\rho }_{min}$ to 1 × 10⁻³ helps to avoid matrix singularity during layer-by-layer strength analysis and static mechanical analysis.

4.2. Sensitivity Analysis

To conduct stress-based topology optimization, it is necessary to provide the sensitivities of the global stress measure on the element densities. The Method of Moving Asymptotes (MMA) [51] is employed to solve the optimization problem, which requires the first-order sensitivity information of the objective function and constraints. The gradient ${\sigma }_{PN}^i$ is calculated using the chain rule as follows:

$$\frac{\partial {\sigma }_{PN}^i}{\partial {\rho }_j}=\frac{\partial {\sigma }_{PN}^i}{\partial \stackrel{=}{\rho }} \frac{\partial \stackrel{=}{\rho }}{\partial \overline{\rho }}\frac{\partial \overline{\rho }}{\partial \tilde{\rho }}\frac{\partial \tilde{\rho }}{\partial {\rho }_j}$$

(21)

In the equation, $\stackrel{=}{\rho }$ represents the physical density field after the three layers of filtering projection. $\frac{\partial \stackrel{=}{\rho }}{\partial \overline{\rho }}$, $\frac{\partial \overline{\rho }}{\partial \tilde{\rho }}$, and $\frac{\partial \tilde{\rho }}{\partial {\rho }_j}$, respectively, denote the standard modifications of sensitivity due to the AM filter, Heaviside projection filter, and smoothing filter.

For $\frac{\partial {\sigma }_{PN}^i}{\partial \stackrel{=}{\rho }}$, we can obtain

$$\begin{gathered} \frac{\partial {\sigma }_{PN}^i}{\partial \stackrel{=}{\rho }}={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_e}\right)}^T\frac{\partial \eta \left({\stackrel{=}{\rho }}_e\right){\mathit{\sigma }}_e}{\partial \stackrel{=}{\rho }}\right]\right) \\ ={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_e}\right)}^T\frac{\partial \eta \left({\stackrel{=}{\rho }}_e\right)}{\partial \stackrel{=}{\rho }}{\mathit{\sigma }}_{\mathit{e}}\right]\right)+{{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\eta \left({\stackrel{=}{\rho }}_e\right)\frac{\partial {\mathit{\sigma }}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}\right]\right) \end{gathered}$$

(22)

The above equation can be written in the following form:

$$\left\{\begin{array}{c}\frac{\partial {\sigma }_{PN}^i}{\partial {\rho }_j}=T_1+T_2\\ T_1={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\frac{\partial \eta \left({\stackrel{=}{\rho }}_e\right)}{\partial \stackrel{=}{\rho }}{\mathit{\sigma }}_{\mathit{e}}\right]\right)\\ T_2={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\eta \left({\stackrel{=}{\rho }}_e\right)\frac{\partial {\mathit{\sigma }}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}\right]\right)\end{array}\right.$$

(23)

For $\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}$, it can be expressed as follows:

$$\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}={{\displaystyle \sum }}_{e=1}^{nele}\left({({\left(\partial {\widehat{\sigma }}_{vm,e}\right)}^p)}^{\frac{1}{p}-1}{\left(\partial {\widehat{\sigma }}_{vm,e}\right)}^{p-1}\right)$$

(24)

The derivative of the element von Mises stress with respect to its stress components can be expressed as

$$\left\{\begin{array}{c}\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\sigma }}_{ex}}=\frac{1}{2{\widehat{\sigma }}_{vm,e}}\left(2{\widehat{\sigma }}_{ex}-{\widehat{\sigma }}_{ey}\right)\\ \frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\sigma }}_{ey}}=\frac{1}{2{\widehat{\sigma }}_{vm,e}}\left(2{\widehat{\sigma }}_{ey}-{\widehat{\sigma }}_{ex}\right)\\ \frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\tau }}_{exy}}=\frac{3}{{\widehat{\sigma }}_{vm,e}}{\widehat{\tau }}_{exy}\end{array}\right.$$

(25)

The $\frac{\partial \eta \left({\stackrel{=}{\rho }}_e\right)}{\partial \stackrel{=}{\rho }}$ is simply obtained as

$$\frac{\partial \eta \left({\stackrel{=}{\rho }}_e\right)}{\partial \stackrel{=}{\rho }}=q{\stackrel{=}{\rho }}_e{}^{q-1}\frac{\partial {\stackrel{=}{\rho }}_e}{\partial \stackrel{=}{\rho }}$$

(26)

The analytic form of $\frac{\partial {\sigma }_e}{\partial \stackrel{=}{\rho }}$ can be expressed as

$$\frac{\partial {\mathit{\sigma }}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}={\theta }_e^i{\mathit{D}}_0{\mathit{B}}_{\mathit{e}}\frac{\partial {\mathit{u}}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}={\theta }_e^i{\mathit{D}}_0{\mathit{B}}_{\mathit{e}}{\mathit{L}}_{\mathit{e}}\frac{\partial {\mathit{U}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}$$

(27)

Then, we obtain the expression of $T_2$ through substitutions

$$T_2={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}\left[{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\eta \left({\stackrel{=}{\rho }}_e\right)\frac{\partial {\mathit{\sigma }}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}\right]\right)$$

$$={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\eta \left({\stackrel{=}{\rho }}_e\right){\theta }_e^i{\mathit{D}}_0{\mathit{B}}_{\mathit{e}}{\mathit{L}}_{\mathit{e}}\frac{\partial {\mathit{U}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}\right)$$

(28)

Given the unknown $\frac{\partial U_i}{\partial \stackrel{=}{\rho }}$, we take the derivatives on the governing equation and derive

$$\frac{\partial {\mathit{K}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}{\mathit{U}}_{\mathit{i}}+{\mathit{K}}_{\mathit{i}}\frac{\partial {\mathit{U}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}=\frac{\partial {\mathit{F}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}$$

(29)

We solve the above formula and then we have

$$\frac{\partial {\mathit{U}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}={\mathit{K}}_{\mathit{i}}{}^{-1}\left(\frac{\partial {\mathit{F}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}-\frac{\partial {\mathit{K}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}{\mathit{U}}_{\mathit{i}}\right)$$

(30)

Thus, $T_2$ can be expressed as

$$T_2={{\displaystyle \sum }}_{e=1}^{nele}\left(\frac{\partial {\sigma }_{PN}^i}{\partial {\widehat{\sigma }}_{vm,e}}{\left(\frac{\partial {\widehat{\sigma }}_{vm,e}}{\partial {\widehat{\mathit{\sigma }}}_{\mathit{e}}}\right)}^T\eta \left({\stackrel{=}{\rho }}_e\right){\theta }_e^i{\mathit{D}}_0{\mathit{B}}_{\mathit{e}}{\mathit{L}}_{\mathit{e}}{\mathit{K}}_{\mathit{i}}{}^{-1}\left(\frac{\partial {\mathit{F}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}-\frac{\partial {\mathit{K}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}{\mathit{U}}_{\mathit{i}}\right)\right)$$

(31)

where $\frac{\partial F_i}{\partial \stackrel{=}{\rho }}$ is obtained from the previous calculation,

$$\frac{\partial {\mathit{F}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}={\displaystyle \sum }_{e=1}^{nele}\frac{3}{4}{\delta }_e^i{\stackrel{=}{\rho }}_e^2\frac{\partial {\stackrel{=}{\rho }}_e}{\partial \stackrel{=}{\rho }}g{a}^{2}b{\mathit{L}}_{\mathit{e}}^{\mathit{T}}$$

(32)

For $\frac{\partial {\mathit{K}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}$, it can be expressed as

$$\frac{\partial {\mathit{K}}_{\mathit{i}}}{\partial \stackrel{=}{\rho }}={{\displaystyle \sum }}_{e=1}^{nele}{\theta }_e^i{\mathit{L}}_{\mathit{e}}^{\mathit{T}}\frac{\partial {\mathit{K}}_{\mathit{e}}}{\partial \stackrel{=}{\rho }}{\mathit{L}}_{\mathit{e}}={{\displaystyle \sum }}_{e=1}^{nele}{\theta }_e^i{\mathit{L}}_{\mathit{e}}^{\mathit{T}}\left(\frac{\partial {\eta }_e^E}{\partial {\stackrel{=}{\rho }}_e}\frac{\partial {\stackrel{=}{\rho }}_e}{\partial \stackrel{=}{\rho }}{\mathit{K}}_0\right){\mathit{L}}_{\mathit{e}}$$

(33)

In the formula, ${\eta }_e^E$ is the elastic modulus interpolation function, and ${\mathit{K}}_0$ is the element stiffness matrix in full solid state. The detailed expression of the element stiffness is given by the following formula:

$${\mathit{K}}_0=\int {\mathit{B}}^T{\mathit{D}}_0\mathit{B}d{\Omega }_e$$

(34)

5. Numerical Examples

In this section, we employ two classic two-dimensional cases, the MBB beam and cantilever beam, to validate the efficacy of the proposed method. In all numerical implementations, the “45-degree rule” is adopted, i.e., the maximum self-supporting overhang angle is set as 45 degrees. Finite element analysis is performed using four-node quadrilateral elements. The MMA optimizer is utilized as the solver with a move limit of 0.1. The sharpness parameter for the Heaviside projection filter starts from 1, doubles for every 50 iterations, and ends with a maximum value of 32.

Regarding the implementation specifications, finite element analysis is performed on the general minimum compliance topology optimization result with only the volume constraint and no stress constraint. The stress fields due to adhesion loads are simulated to disclose the prone-to-failure features and provide references for setting the stress restriction upper bounds. Thereafter, the stress-constrained and volume-constrained minimum compliance topology optimization is carried out to design the process-safe optimal structures. As stated earlier, setting $q$ to 0.5 can help to avoid gray elements remaining in the optimized structure for the traditional stress-constrained SIMP setup. However, artificial high stresses are likely to appear at the boundary non-solid elements since the AM filter is applied after Heaviside projection and the approximated $\mathrm{smax}$ and $\mathrm{smin}$ operators unavoidably lead to non-solid elements. $q=0.5$ cannot eliminate these non-solid effects; instead, it artificially magnifies the local stresses that disturb the convergence. Hence, in the following case study, $q$ is set to 2 to prevent non-physical stress concentration.

5.1. The MBB Beam Structure

The first case is an MBB beam example, as illustrated in Figure 4. Due to the symmetry condition, only the right half of the structure is considered in optimization with the discretization of $nelx=180$ and $nely=60$. A concentrated vertical force ($F=1$) is loaded at the top left corner of the half domain, while the bottom right corner is supported on a roller. The left edge adopts the symmetry boundary condition. The element size of 1 × 1 is adopted. The volume fraction upper limit is chosen to be 0.5. The filter radius $rmin$ is set to 5. Other parameters are selected as follows: the SIMP interpolation parameter $p\_SIMP$ is 5; the RAMP interpolation parameter $p\_RAMP$ is 20; the P-norm parameter $p$ is 10.

Using SIMP interpolation, the minimum compliance of 194.8950 is first obtained without applying the stress constraint, and if we apply the layer-based adhesion load to this result, the maximum P-norm stress occurs at the 49th layer’s adhesion load, reaching 14.9651. Hence, to reduce the risk of failure, we take the stress limits of 10, 5, and 3, respectively, to perform the process of stress-constrained topology optimization, obtaining the optimized structures and corresponding stress profiles for the most stress-concentrated load steps, as presented in Figure 5. Note that the full MBB structure is demonstrated instead of only a symmetric half.

With the P-norm stress limit of 10, the optimized structural compliance is 194.4829 and the maximum P-norm stress occurs at the 49th layer’s load step, reduced to 8.4181.

With the P-norm stress limit of 5, the optimized structural compliance is 208.4969 and the maximum P-norm stress occurs at the 37th layer’s load step, reduced to 4.2599.

With the P-norm stress limit of 3, the optimized structural compliance is 222.9046 and the maximum P-norm stress occurs at the 53rd layer’s adhesion load, reduced to 3.0000.

The convergence curves under the three P-norm stress constraint conditions are shown in Figure 6.

From the iteration curves in Figure 6, we can see that the convergence of the iteration curves is difficult with the SIMP interpolation. The curves fluctuate violently and the fluctuation intensity increases with the enhanced strength requirement.

Using RAMP interpolation, the minimum compliance of 187.7462 is obtained without the stress constraint, and if we apply the layer-based adhesion load, the maximum P-norm stress occurs at the 47th layer’s adhesion load, reaching 18.3128. Hence, to reduce the failure risk, we again take the P-norm stress limits of 10, 5, and 3, respectively, to perform the process of stress-constrained topology optimization, obtaining the optimized structures and corresponding stress profile plots for the most stress-concentrated load steps, as shown in Figure 7.

With the stress limit of 10, the optimized structural compliance is 186.1907 and the maximum P-norm stress occurs at the 47th layer’s adhesion load with the peak value of 8.7847.

With the stress limit of 5, the optimized structural compliance is 187.4594 and the maximum P-norm stress occurs at the 50th layer’s adhesion load with the peak value of 4.9952.

With the stress limit of 3, the optimized structural compliance is 190.3872 and the maximum P-norm stress occurs at the 50th layer’s adhesion load with the peak value of 2.9935.

Figure 8 summarizes the optimized structural compliances subject to different constraint conditions and interpolation methods. It can be seen that under the currently selected parameters, using RAMP interpolation leads to structures with better overall stiffness compared to SIMP interpolation. However, the difference between the two interpolation methods may be decreased by adjusting the interpolation details.

The convergence curves under the three stress constraint conditions for the RAMP interpolation scenario are shown in Figure 9.

From Figure 9, we can see that the iteration curves converge faster with the RAMP interpolation, and fewer fluctuations are observed than in Figure 6. Additionally, the overall stiffness deterioration due to involving stress constraints is moderate for the RAMP interpolation.

5.2. Cantilever Beam

The second case is performed on a cantilever beam, as illustrated in Figure 10, with $nelx=160$ and $nely=80$. The load is applied on the bottom right corner, while the left edge is completely fixed. The volume fraction upper limit is chosen to be 0.5. The filter radius $r_{min}$ is set to 5. Other parameters are selected as follows: the SIMP interpolation parameter $p\_SIMP$ is 5; the RAMP interpolation parameter $p\_RAMP$ is 20; the P-norm parameter $p$ is 10.

(1)

Using SIMP interpolation, the minimum compliance of 81.5803 is first obtained without imposing the stress constraint, and if we apply the layer-based adhesion load to this result, the maximum P-norm stress occurs at the 59th layer’s adhesion load, reaching 17.906. Hence, to reduce the risk of failure, we take the stress limits of 10, 5, and 3, respectively, to perform the process of stress-constrained topology optimization, obtaining the optimized structures and corresponding stress profiles for the most stress-concentrated load steps, referring to Figure 11.

With the P-norm stress limit of 10, the optimized structural compliance is 76.3174 and the maximum P-norm stress of 10.0000 occurs at the 45th layer’s adhesion load step.

With the P-norm stress limit of 5, the optimized structural compliance is 76.5998 and the maximum P-norm stress of 4.9990 occurs at the 60th layer’s adhesion load step.

With the P-norm stress limit of 3, the optimized structural compliance is 79.3136 and the maximum P-norm stress of 3.0000 occurs at the 47th layer’s adhesion load step.

The convergence curves under the three stress constraint conditions are given in Figure 12.

From the iteration curves in Figure 12, we can see that when using SIMP interpolation, compared to the MBB beam case, the iteration curves are smoother and converge more easily. The stiffness of the optimized results gradually deteriorates with the increased strength requirement, but the overall stiffness change is small for the stress limits of 10 and 5. It is unexpected that the optimization result without applying the stress constraint has the lowest stiffness, possibly due to the local optimum issue.

(2)

Using RAMP interpolation, the minimum compliance of 68.2931 is obtained without the stress constraint, and if we apply the layer-based adhesion load, the maximum stress occurs at the 63rd layer’s adhesion load step, reaching 17.1274. Hence, to reduce the failure risk, we again take the stress limits of 10, 5, and 3, respectively, to perform the process of stress-constrained topology optimization, obtaining the optimized structures and corresponding stress profiles for the most stress-concentrated load steps, shown in Figure 13.

With the P-norm stress limit of 10, the optimized structural compliance is 68.6842 and the maximum P-norm stress of 9.8595 occurs at the 62nd layer’s adhesion load.

With the P-norm stress limit of 5, the optimized structural compliance is 68.7421 and the maximum P-norm stress of 4.9990 occurs at the 64th layer’s adhesion load.

With the P-norm stress limit of 3, the optimized structural compliance is 69.7818 and the maximum P-norm stress of 2.9969 occurs at the 65th layer’s adhesion load.

Figure 14 summarizes the optimized structural compliances subject to different constraint conditions and interpolation methods. The observation remains the same as the MBB beam case, namely that under the current parameter selection, the RAMP interpolation produces structures with better global stiffness.

The convergence curves under the three stress constraint conditions are given in Figure 15. Similar to the first numerical example, the stiffness condition of the structure gradually deteriorates with the increased strength requirement, and the overall stiffness deterioration is evidently milder than with the SIMP interpolation.

6. Experimental Validation

In order to validate the efficacy of the proposed method, an experimental test is designed and carried out in this section. An MBB beam with an 8:1 aspect ratio is optimized and manufactured for this experimental investigation, since the increased structure details provide a higher possibility of failure. The stress-constrained topology optimization is performed with SIMP interpolation and the following parameters: $nelx=240$, $nely=60$, 50% volume fraction upper limit, SIMP penalization factor $p\_SIMP$=4, filter radius of 4, P-norm parameter $p=10$, stress interpolation factor $q=2$, 60 layers of adhesion load, and the MMA optimizer move limit of 0.1.

Finite element analysis on the topology optimization result without the stress constraint shows the maximum P-norm stress occurring at the 52nd layer’s adhesion load, reaching 38.01. The overall compliance of the structure is 399.5586. The 3D model for additive manufacturing obtained after post-processing the optimization result is shown in Figure 16.

Then, we set a process stress upper limit of 2 and obtain the optimized design as shown in Figure 17. The maximum stress occurs at the 34th layer’s adhesion load, having the peak P-norm stress of 2.000. The overall structural compliance is 412.1184, slightly increased compared to the non-constrained design in Figure 15. However, the weak and stress-concentrated thin struts are successfully removed in this new design.

As shown in Figure 18, a thickness of 0.5 mm is selected for model generation. Three repeated of print attempts are made for each of the above two compared designs. The print results are demonstrated in Figure 19. It is observed that the regular topology optimization designs all fail at the stress-concentrated thin strut, while the process stress-constrained designs are all successfully printed with complete structural details. Therefore, the proposed method enhances the in-process structural strength and can resolve the issue of print failure due to the repeating adhesion loads.

7. Conclusions

This paper has developed a topology optimization design method for DLP additive manufacturing that addresses the process-related structural strength issue, leading to the following conclusions.

(1)

To mimic the in-process stress concentrations, an SIMP-like interpolation is proposed to simulate the adhesion forces between the workpiece and resin base for DLP additive manufacturing.

(2)

A greater than 1 penalization parameter is adopted for the stress term to prevent artificial high stresses at the boundary non-solid elements since the AM filter is applied after Heaviside projection and the approximated $\mathrm{smax}$ and $\mathrm{smin}$ operators unavoidably lead to non-solid elements.

(3)

Both the SIMP and RAMP interpolations allow us to achieve topology optimization designs considering the maximum volume fraction and P-norm stress constraints. SIMP interpolation exhibits more fluctuations in converging the optimization process. RAMP interpolation has a smoother convergence process, although some cases show a small amount of irremovable gray elements.

(4)

Experiments validate the necessity and effect of imposing the stress constraints. The stress-concentrated thin struts are replaced by thickened structural features through size and shape enhancement, and the overall stress field tends to distribute more evenly. Hence, failures are not encountered in the modified designs and, simultaneously, the local shaping accuracy is guaranteed as well.

(5)

The layerwise P-norm stresses are summarized in Figure 20 for both the SIMP and RAMP interpolations. Apparently, a stricter stress restriction leads to an overall reduction in P-norm stresses, indicating the effectiveness of the layerwise stress constraints.

(6)

A stricter stress restriction, i.e., lower P-norm stress upper limit, causes the decreased overall stiffness of the optimized structure for both the SIMP and RAMP interpolations.

(7)

Interestingly, topological optimization that considers the maximum volume fraction and stress constraints may yield better structural stiffness than topological optimization that only considers the volume constraints. This phenomenon deserves careful further investigation.

(8)

The build orientation plays an equally important role in topology optimization design considering the process stress constraints in additive manufacturing. Appropriate build directions may allow us to resolve the maximum P-norm stress constraints without significantly impacting the structural stiffness. This topic will be investigated in our forthcoming research.

(9)

Finally, considering the computational cost, the finite element analysis simulating the layer-based adhesion loading conditions and the stress term-related adjoint equations occupies the majority of the calculation time. Adopting parallel computing could greatly reduce the computing time, especially for larger-scale problems, and this scheme deserves further exploration as well.