Deep energy method in topology optimization applications

Deep neural networks (DNNs) consist of multiple layers of interconnected neurons, and the layers themselves are also connected. The configuration and function of DNNs bear resemblance to that of a human brain. Fig. 1 shows a typical structure of a fully connected neural network model.

The architecture of a DNN can be characterized by the number of hidden layers and the number of neurons in each hidden layer. The neurons of consecutive layers are connected by a set of weights $\bm{W}$ and biases $\bm{b}$. The output $\bm{y}^{i}$ of layer $i$ can be related to its corresponding weights and biases as:

where $f_{\rm{act}}^{i}$ denote the activation function for this layer. Training of the DNN refers to the iterative process where the weights and biases of the model are updated by means of gradient descent [39] such that the value of a user-defined loss function is minimized.

In this section, we describe DEM in the context of the elasticity equation in structural mechanics. We consider a homogeneous, isotropic, linear elastic body undergoing small deformation. We remark that DEM is not limited to linear elastic materials and has been applied in nonlinear hyperelastic material models [27]. In the absence of any body and inertial forces, the equilibrium equation in tensorial form reads:

The system is subject to Dirichlet boundary condition:

and/or Neumann boundary condition:

In the small deformation setting, the strain tensor is given by:

Stress can be computed from the constitutive law as:

where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio. DEM seeks the solution to the equilibrium equation via the principle of minimum potential energy. For a body in static equilibrium with no applied body forces, the potential energy of the system reads:

The loss function ($\mathcal{L}$) in DEM is defined identically as the potential of the system:

The training of the DEM model $\mathcal{M}:\mathbb{R}^{n}\to\mathbb{R}^{n}/\bm{X}\Rightarrow\tilde{\bm{u}}(\bm{X})$ ¹¹1Note that the immediate output of DEM model is $\tilde{\bm{u}}$, which is different from $\bm{u}$ in Eq. (5) and Eq. (7). can be viewed as minimizing the loss function (system potential energy) $\mathcal{L}$ through an optimizer. The solution to the elasticity problem, as given by the DEM model, is defined as:

In this work, the L-BFGS optimizer [40] was used to train the model. The learning rate is another parameter that can be changed by the user to achieve optimal framework performance. Similar to the early stopping schedule commonly used in model training, we monitor the relative change of the loss function and stop the training once the relative change in loss function value is less than the user-specified tolerance $\epsilon_{tol}$.

The underlying architecture of the DEM model is a fully connected feed-forward network, which is similar to the one depicted in Fig. 1 and was described in the work of Nguyen et al. [27]. Fourier transform was applied using the package Random Fourier Features Pytorch (RFF) [41] to the input features to transform them to frequency domain, which was shown to increase the accuracy of PINNs in the work of Wang et al. [42]. The NN architecture is fully parametrizable, and is characterized by the number of layers, number of neurons (per layer), the activation function, and two standard deviation values used to initialize the weights and biases of the MLP and RFF. Hyperparameter optimization can be performed to improve the performance of the DEM model.

In the context of elasticity problems, three different types of boundary conditions (BCs) are common, namely the Neumann, Dirichlet, and periodic BCs. Enforcement of Neumann BCs is straightforward. The value of the prescribed traction $\bar{\bm{t}}$ enters directly in the boundary integral part of Eq. (7), while the boundary displacements are provided by the DEM model. Traction-free surfaces are satisfied trivially as the boundary integral vanishes on these surfaces. A surface (in 3D) or edge (in 2D) integration rule needs to be defined for numerically evaluating the surface integral in regions where traction is prescribed.

Dirichlet BCs are typically enforced via the penalty method by adding a penalty term to the loss function [27, 43, 17]:

leading to a modified loss function definition:

where $w$ is a user-defined weight parameter. To avoid defining the additional parameter $w$, we adopted a direct approach to enforce Dirichlet BCs via an additive decomposition:

where $*$ denotes element-wise multiplication between two vectors, the vectors $\bm{m}$ and $\bm{u}_{0}$ are chosen such that:

By design, $\bm{u}(\bm{X})$ satisfies all Dirichlet boundary conditions for all arbitrary $\tilde{\bm{u}}$ produced by the DEM model.

Periodic BCs are commonly used when simulating a representative volume element that repeats infinitely in space. It is also used when simulating unit cells of periodic meta materials. In the context of the FEM, periodic BCs are typically stated in the form of constraint equations [44]:

where $k-$ and $k+$ denote the left and right instances of the $k^{th}$ periodic boundary pair, and $\epsilon^{0}_{ij}$ is the (known) applied strain tensor. For easier convergence of the DEM training process, we adopted a mixed implementation for periodic BCs, where we treat zero and non-zero applied strain components differently. The boundary displacement induced by a non-zero applied strain component is enforced as a pair of Dirichlet BCs, and zero applied strain components are enforced by penalty. Without loss of generality, let the non-zero component of $\epsilon^{0}_{ij}$ be $\epsilon_{a}$, which is applied on the $k^{th}$ periodic boundary pair. This leads to two Dirichlet boundary conditions:

where $\Delta L_{k}$ is the distance between the $k+$ and $k-$ boundary pair. For all other boundary pairs that correspond to zero applied strain components, the constraint equation Eq. (14) reduces to:

Therefore, the following penalty term is added to the loss function:

where $E$ is the Young’s modulus of the material (a physics-informed penalty weight) and $N_{p}$ is the number of points of the $k^{th}$ periodic boundary.

In this section, we describe topology optimization (TO) in the context of DEM. The mathematical formulation of the TO problem using DEM can be stated as:

where $f(\bm{u},\bm{\rho})$, $\bm{\rho}$, $V_{e}$ and $\bar{V}$ denote the objective function to minimize, element density vector, element volume and target volume, respectively. A typical objective of TO is to minimize the system compliance subject to a volume/mass constraint. For compliance minimization, the objective function is:

where $\bm{\sigma}^{*}(\bm{\rho})$ is the stress computed from the classical solid isotropic material penalization (SIMP) scheme [45]:

The penalty exponent $p$ is taken to be 3 in this work.

For designing periodic 2D meta structures, a commonly employed objective is to maximize the homogenized shear modulus of the unit cell [31, 44, 46]:

where $A$ is the area of the 2D element, $\bm{\sigma}^{p}$ and $\bm{\epsilon}^{p}$ denote the stress and strain fields computed from the periodic displacement field subjected to the following applied simple shear loading ($\bm{\epsilon}^{0}$ in Eq. (14) ):

To fit into a minimization framework, the objective function is defined as:

Due to the introduction of the element density vector, we correspondingly modify the definition of the potential energy and loss function as:

In this work, we chose to evaluate the volume integral in Eq. (24) by Gauss quadrature integration. This necessitates the formation of quadrilateral (in 2D) or hexahedral (in 3D) isoparametric finite elements from the structured nodes $\bm{X}$ in the DEM domain, which is very similar to a FEM-type treatment. The need to form elements is a major disadvantage of the framework, as compared to other meshless DEM methods like those developed in [17, 27] which only require nodes but not elements to exist in the simulation domain. The gradient operator in Eq. (5) is evaluated by the gradients of the finite element shape functions, and not from automatic differentiation of the DEM model. The density vector $\bm{\rho}$ is defined at the cell center, identical in classical FEM-based TO. The method of moving asymptotes (MMA) [47] is used to iteratively update the density based on the sensitivity.

In this section, we describe the sensitivity analysis in the context of DEM. Comparing Eq. (19) and Eq. (23), we see that both are similar in form to a classical compliance minimization problem in FEM-based TO. Since the compliance minimization problem is self-adjoint, the adjoint vector is simply the negative of the displacement, and therefore, no additional adjoint analysis is needed to evaluate the sensitivity. In this case, one simply needs to express the sensitivity in terms of the displacement field and its gradients produced by DEM. For compliance minimization problem (Eq. (19)), the sensitivity of element $e$ is given by:

Similarly, for maximization of shear modulus in 2D meta material unit cells (Eq. (23)), the sensitivity of element $e$ is given by:

To avoid checkerboarding and mesh-dependence, the density filtering technique introduced by Bruns et al. [48] was employed. In this approach, the MMA solver does not directly manipulate $\bm{\rho}$, but instead manipulates on a pseudo-density (the design variables) $\bm{\xi}$ that is related to $\bm{\rho}$ as:

where $\bm{W}$ is the normalized filter matrix whose rows are given by the normalized weights $\bar{q}_{ij}$:

Algorithm 1 summarizes the process of TO using the DEM method.

Two different cases are presented in this example. In the first case, a 2D rectangular bridge of size 12-by-2 m² was subjected to a downward surface traction at the center of the top edge over a length of 0.5 m. The bridge was fixed at its left and right edges and 3751 nodes were placed into the domain forming a 121-by-31 grid. In the second case, a 2D rectangular beam of size 10-by-5 m² was subjected to a downward concentrated load at the center of the right edge, while the left edge was fixed. 4186 nodes were placed into the domain, forming a 91-by-46 grid. The relative convergence tolerance $\epsilon_{tol}$ was set to $5\times 10^{-6}$ for all DEM training. As comparison, FEM-based TO was applied to solve the same problem using identical node layout.

The designs from DEM- and FEM-based TO at different design iterations are compared in Figures 2 and 3. The relative compliance reduction for both cases as the TO optimizer progressed is presented in Fig. 4a. To give a quantitative comparison of the similarity between DEM- and FEM-based designs, the dice similarity coefficient (DSC) was computed for both cases and is presented in Fig. 4b. The DSC is defined as:

where $\bm{I}$ represents the binarized density using a (dynamic) threshold of 0.5 times the current maximum density. Finally, the computational cost for running the two cases using DEM and FEM²²2The highly optimized Python TO code found in [52] was used for time comparison, which was based on the MATLAB code by Andreassen et al. [53] are summarized in Fig. 4c and are compared in Table 2.

From Fig. 2 and Fig. 3, we clearly see that the DEM- and FEM-based TO frameworks produced almost identical simulation design evolution history and final designs. The similarity in designs is further demonstrated in Fig. 4b, where we see that the DSC is greater than 90% for the final designs in both cases. For the bridge case, we noticed that the DSC remains greater than 98% for all design iterations, showing close resemblance between the DEM- and FEM-based approaches. The two frameworks also produce designs of similar performance, as is evident in Fig. 4a, where we see that the compliance reduction history is almost identical between the two methods due to the use of identical MMA optimizer. Nonetheless, Table 2 shows that FEM still holds a massive advantage over DEM in terms of computational time, which agrees with the trends observed in the similar DEM-based TO work like [38]. However, a very tight relative tolerance of $5\times 10^{-6}$ was used in training, which could be loosened to gain computational efficiency. From Fig. 4c, we highlight two important features of DEM-based TO. For both cases, we noticed that the training time shows a decreasing trend as design progressed. This is due to the fact that the weights and biases from the previous TO iteration are used as the initial condition for the next iteration, a form of transfer learning. This leads to gradual training time reduction. This is in sharp contrast of FEM-based TO, where the computational time for each iteration is constant, and can be very beneficial when many TO iterations are needed. Also, we note that the DEM training time went up by merely 5.4% from case 1 to case 2, where the FEM-based TO time went up over 200%. The insensitivity of training time to domain size is due to the fact that the number of parameters of DEM-based TO is controlled by the DNN architecture, not by the number of nodes used to discretize the domain, which is beneficial for problems that need to be resolved by a fine mesh.

The second example concerns a 3D rectangular bridge of size 12-by-2-by-2 m³. It was fixed at both end surfaces and was subjected to a downward load at the center of the top surface in a square region of size 0.5-by-0.5 m². 75625 nodes were placed into the domain, forming a 121-by-25-by-25 grid. $\epsilon_{tol}$ was set to $5\times 10^{-5}$. As comparison, Abaqus/Standard [54] was used to perform TO using the built-in TO functionalities and optimizer.

The final designs from DEM and FEM are presented in Fig. 5a and Fig. 5b, respectively. The STL file for the final FEM design was extracted using a density threshold of 0.8, and the same threshold was used to obtain the voxelated DEM final design. The two designs are overlaid and presented in Fig. 5c and Fig. 5d. The relative compliance reduction and training time was shown in Fig. 6. The total TO time for DEM was 1050.4s.

From the overlaid plot and the section view in Fig. 5, we see that the DEM-based design largely resembles the FEM-based design, especially near both ends of the bridge. Some design differences can be seen in the center, especially from the section view. From Fig. 6a, we see that DEM was able to quickly reduce the compliance in less than 10 iterations, but has a final compliance that is larger than that of the FEM design. The evolution histories are different, which roots from the fact that different optimizers were used to optimize the material distribution. From Fig. 6b, we again see that training time decreased as design progressed, a trend that is consistent with Section 3.1. Finally, we highlight that the same hyperparameters in Section 3.1 were used directly in this 3D example without any additional hyperparameter optimization, which shows the robustness of the DEM model with respect to hyperparameter selection.

In this example, a 2D unit cell of size 10-by-10 m² was subjected to a shear load at the top edge. The unit cell was subjected to periodic BCs in the X- and Y-direction. Further, 6561 nodes were placed into the domain forming a 81-by-81 grid. The objective is to maximize the homogenized shear modulus while maintaining the volume fraction at 40%. It is well-known that the final optimized design depends heavily on the initial material distribution [31, 46, 44]. Therefore, three different initial designs were considered with different initial hole diameters, as shown in the first row of Fig. 7. $\epsilon_{tol}$ was set to $5\times 10^{-6}$. For comparison, the DEM-based TO results are compared to their FEM-based counterparts produced by the MATLAB [55] code developed by Xia et al. [44].

The initial and final designs from DEM- and FEM-based TO are compared in Fig. 7. The relative change in shear modulus for both cases over all the iterations are presented in Fig. 8a. Further, DSC coefficient is calculated for the binarized final designs obtained from both cases using Eq. (29) and is presented in Fig. 8b. The DEM training time for the three cases were 579.7, 934.4 and 1125.7s, respectively. Training time is summarized in Fig. 8c.

From Fig. 7, we clearly see that the DEM- and FEM-based TO frameworks produced very similar final designs for all three initial designs. However, comparing Fig. 7f and Fig. 7i, we noticed an important deficiency of the DEM-based TO method. The final FEM-based TO design (Fig. 7i) has four-fold symmetry, while the DEM-based design (Fig. 7f) does not. In fact, the DEM-based final design is not periodic at the lower left corner. This difference is due to that fact that periodic BCs are enforced in FEM-based simulations by constraining the corresponding degrees of freedom of the global stiffness matrix, hence periodicity of the displacement field (hence the material distribution) is also guaranteed. Also, if the initial design is symmetric, symmetry persists in subsequent designs. This is not the case for DEM-based TO, as the periodic BCs were enforced by penalty. The periodic penalty loss gradually creep up as the design progressed, and the designs gradually lost periodicity due to minor violations of periodicity in the displacement field and material distribution. Nonetheless, it is observed that the DSC is greater than 90% for the final designs in all three cases. The low DSC during the initial stage of evolution is a result of the different optimizers used in the DEM- and FEM-based TO code. Similar to the two examples presented above, we again see gradual reduction of training time as TO iterations progressed.