Topology optimization for transient response of structures subjected to dynamic loads

TL;DR

This paper develops a topology optimization framework for structures under transient dynamic loads, enabling tailored responses through advanced numerical methods and material interpolation techniques.

Contribution

It introduces a gradient-based optimization approach using MMA, adjoint sensitivities, and stabilization methods for transient structural problems.

Findings

Optimized designs effectively tailor transient responses.

The framework handles both static and dynamic loading cases.

Numerical stabilization improves manufacturability of designs.

Abstract

This paper presents a topology optimization framework for structural problems subjected to transient loading. The mechanical model assumes a linear elastic isotropic material, infinitesimal strains, and a dynamic response. The optimization problem is solved using the gradient-based optimizer Method of Moving Asymptotes (MMA) with time-dependent sensitivities provided via the adjoint method. The stiffness of materials is interpolated using the Solid Isotropic Material with Penalization (SIMP) method and the Heaviside Projection Method (HPM) is used to stabilize the problem numerically and improve the manufacturability of the topology-optimized designs. Both static and dynamic optimization examples are considered here. The resulting optimized designs demonstrate the ability of topology optimization to tailor the transient response of structures.

Tables

Table 1. Table 1: Model and material parameters for the cantilever beam model.

Description	Parameter	Value	Units
minimum allowable radius	$r_{min}$	$1.375\times {10}^{-2}$	m
elemental length	$h_l$	$3.333\times {10}^{-3}$	m
minimum allowable density	${\rho }_{min}$	$1.0\times {10}^{-4}$
Young’s modulus of pure solid	$E_0$	$2.0\times {10}^{11}$	Pa
Poisson’s ratio	$\nu$	0.33
beam length	$L$	0.9	m
beam thickness	$t_b$	0.01	m
maximum applied force	$f_{\max }$	$1.0\times {10}^4$	N
maximum allowable time	$t_f$	$\pi$	s
time step size	$\Delta t$	$\pi /20$	s
Heaviside projection parameter	$\beta$	50
SIMP exponential penalty term	$\eta$	1,…,10
lower bound of the design variables	${\varphi }_{min}$	0.0
upper bound of the design variables	${\varphi }_{max}$	1.0
maximum allowable volume fraction	$V_{max}$	0.35

Equations

$$\begin{array}{r}\mathbf{M}{\ddot{\mathbf{u}}}^n+\mathbf{C}{\dot{\mathbf{u}}}^n+\mathbf{K}{\mathbf{u}}^n-{\bar{\mathbf{f}}}^n=0,\end{array}$$

$$\begin{array}{r}{\dot{\mathbf{u}}}^n=\left(1-\frac{\tilde{\gamma }}{\tilde{\beta }}\right){\dot{\mathbf{u}}}^{n-1}+\left(1-\frac{\tilde{\gamma }}{2\tilde{\beta }}\right)\Delta t{\ddot{\mathbf{u}}}^{n-1}+\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}\left({\mathbf{u}}^n-{\mathbf{u}}^{n-1}\right),\end{array}$$

$$\begin{array}{r}{\ddot{\mathbf{u}}}^n=-\frac{1}{\tilde{\beta }\Delta t}{\dot{\mathbf{u}}}^{n-1}+\left(\frac{1}{2\tilde{\beta }}-1\right){\ddot{\mathbf{u}}}^{n-1}+\frac{1}{\tilde{\beta }\Delta t^2}\left({\mathbf{u}}^n-{\mathbf{u}}^{n-1}\right),\end{array}$$

$$\begin{array}{r}{\mathbf{R}}^n=\left[\left(\frac{1}{2\tilde{\beta }}-1\right)\mathbf{M}+\left(1-\frac{\tilde{\gamma }}{2\tilde{\beta }}\right)\Delta t\mathbf{C}\right]{\ddot{\mathbf{u}}}^{n-1}-\left[\frac{1}{\tilde{\beta }\Delta t}\mathbf{M}-\left(1-\frac{\tilde{\gamma }}{\tilde{\beta }}\right)\mathbf{C}\right]{\dot{\mathbf{u}}}^{n-1}\\ +\left[\frac{1}{\tilde{\beta }\Delta t^2}\mathbf{M}+\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}\mathbf{C}\right]\left({\mathbf{u}}^n-{\mathbf{u}}^{n-1}\right)+\mathbf{K}{\mathbf{u}}^n-{\bar{\mathbf{f}}}^n\forall n=1,...,N_t.\end{array}$$

$$\begin{array}{r}{\mathbf{R}}^0={\mathbf{u}}^0-{\mathbf{u}}_{\mathrm{int}},\end{array}$$

$$\begin{array}{r}\underset{\mathbf{\varphi }}{\min }Z(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^0,...,{\mathbf{u}}^{N_t})=\sum _{n=0}^{N_t}{z}^n(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^n)\end{array}$$

$$\text{subject to}\{\begin{array}{l}{H}_i(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^0,...,{\mathbf{u}}^{N_t})=\sum _{n=0}^{N_t}{h}_i^n(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^n)=0\forall i=1,...,N_h\\ G_j(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^0,...,{\mathbf{u}}^{N_t})=\sum _{n=0}^{N_t}{g}_j^n(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^n)\le 0\forall j=1,...,N_g\\ {\mathbf{R}}^n(\mathbf{\rho }(\mathbf{\varphi }),{\mathbf{u}}^n,{\mathbf{u}}^{n-1})=0\\ {\varphi }_{min}\le {\varphi }_k\le {\varphi }_{max}\forall k=1,...,N_{\varphi }\end{array},$$

$$\begin{array}{r}\tilde{\mathbf{c}}=[c^0,...,c^{N_t}{]}^T,\\ \tilde{\mathbf{R}}=[{\mathbf{R}}^0,...,{\mathbf{R}}^{N_t}{]}^T,\\ \tilde{\mathbf{u}}=[{\mathbf{u}}^0,...,{\mathbf{u}}^{N_t}{]}^T,\end{array}$$

$$\begin{array}{r}\frac{d\tilde{\mathbf{c}}}{d{\varphi }_k}=\frac{\partial \tilde{\mathbf{c}}}{\partial {\varphi }_k}+\frac{\partial \tilde{\mathbf{c}}}{\partial \tilde{\mathbf{u}}}\frac{d\tilde{\mathbf{u}}}{d{\varphi }_k}.\end{array}$$

$$\frac{d\tilde{\mathbf{R}}}{d{\varphi }_k}=\frac{\partial \tilde{\mathbf{R}}}{\partial {\varphi }_k}+\frac{\partial \tilde{\mathbf{R}}}{\partial \tilde{\mathbf{u}}}\frac{d\tilde{\mathbf{u}}}{d{\varphi }_k}=0.$$

$$\frac{d\tilde{\mathbf{c}}}{d{\varphi }_k}=\frac{\partial \tilde{\mathbf{c}}}{\partial {\varphi }_k}+{\tilde{\mathbf{\lambda }}}^T\frac{\partial \tilde{\mathbf{R}}}{\partial {\varphi }_k},$$

$${\left(\frac{\partial \tilde{\mathbf{R}}}{\partial \tilde{\mathbf{u}}}\right)}^T\tilde{\mathbf{\lambda }}=-{\left(\frac{\partial \tilde{\mathbf{c}}}{\partial \tilde{\mathbf{u}}}\right)}^T.$$

$$\frac{\partial {\mathbf{R}}^0}{\partial {\mathbf{u}}^j}=\{\begin{array}{ll}\mathbf{I}& \forall j=0\\ 0& \forall j=1,...,N_t\end{array},$$

$$\frac{\partial {\mathbf{R}}^n}{\partial {\mathbf{u}}^j}=\{\begin{array}{ll}-\frac{1}{\tilde{\beta }\Delta t^2}\mathbf{M}-\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}\mathbf{C}& \forall j=n-1\\ \frac{1}{\tilde{\beta }\Delta t^2}\mathbf{M}+\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}\mathbf{C}+\mathbf{K}& \forall j=n\\ 0& \forall j\in \left\{1,...,N_t\right\}\backslash \left\{n-1,n\right\}\end{array}.$$

$$\begin{array}{r}\frac{\partial {\mathbf{R}}^n}{\partial {\varphi }_k}=\left[\left(\frac{1}{2\tilde{\beta }}-1\right)\frac{\partial \mathbf{M}}{\partial {\varphi }_k}+\left(1-\frac{\widetilde{\stackrel{~}{\gamma }}}{2\tilde{\beta }}\right)\Delta t\frac{\partial \mathbf{C}}{\partial {\varphi }_k}\right]{\ddot{\mathbf{u}}}^{n-1}\\ -\left[\frac{1}{\tilde{\beta }\Delta t}\frac{\partial \mathbf{M}}{\partial {\varphi }_k}-\left(1-\frac{\tilde{\gamma }}{\tilde{\beta }}\right)\frac{\partial \mathbf{C}}{\partial {\varphi }_k}\right]{\dot{\mathbf{u}}}^{n-1}\\ +\left[\frac{1}{\tilde{\beta }\Delta t^2}\frac{\partial \mathbf{M}}{\partial {\varphi }_k}+\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}\frac{\partial \mathbf{C}}{\partial {\varphi }_k}\right]\left({\mathbf{u}}^n-{\mathbf{u}}^{n-1}\right)+\frac{\partial \mathbf{K}}{\partial {\varphi }_k}{\mathbf{u}}^n.\end{array}$$

$$\begin{split}E^{e}(\rho)=\big{(}(\rho^{e})^{\eta}+\rho_{min}^{e}\big{)}E_{0}^{e},\end{split}$$

$${\mu }^e(\varphi )=\frac{\sum _{j\in {\Omega }_w}{w}_j{\varphi }_j}{\sum _{j\in {\Omega }_w}{w}_j},\text{with}{w}_j=\{\begin{array}{ll}\frac{r_{min}-\parallel x_j-x^e\parallel }{r_{min}}& \text{ if }x\in {\Omega }_w\\ 0& \text{otherwise}\end{array},$$

$${\rho }^e(\varphi )=1-e^{-\beta {\mu }^e(\varphi )}+\frac{{\mu }^e(\varphi )}{{\varphi }_{max}}{e}^{-\beta {\varphi }_{max}},$$

$$\begin{split}\underset{\boldsymbol{\phi}}{\text{min}}&\quad\quad\displaystyle Z=\boldsymbol{u}^{T}\boldsymbol{K}(\phi)\boldsymbol{u}\\ \text{subject to}&\qquad\displaystyle\sum_{e\in\Omega}\rho^{e}(\phi)v^{e}\leq V_{max}\\ &\qquad\displaystyle\boldsymbol{K}(\phi)\boldsymbol{u}=f_{\mathrm{max}}\Big{|}_{A}+f_{\mathrm{max}}\Big{|}_{B}\\ &\qquad\displaystyle\phi_{min}\leq\boldsymbol{\phi}\leq\phi_{max}\end{split},$$

$$\begin{split}\underset{\boldsymbol{\phi}}{\text{min}}&\quad\quad\displaystyle Z=(\boldsymbol{u}_{1})^{T}\boldsymbol{K}(\phi)\boldsymbol{u}_{1}+(\boldsymbol{u}_{2})^{T}\boldsymbol{K}(\phi)\boldsymbol{u}_{2}\\ \text{subject to}&\qquad\displaystyle\sum_{e\in\Omega}\rho^{e}(\phi)v^{e}\leq V_{max}\\ &\qquad\displaystyle\boldsymbol{K}(\phi)\boldsymbol{u}_{1}=f_{\mathrm{max}}\Big{|}_{A}\\ &\qquad\displaystyle\boldsymbol{K}(\phi)\boldsymbol{u}_{2}=f_{\mathrm{max}}\Big{|}_{B}\\ &\qquad\displaystyle\phi_{min}\leq\boldsymbol{\phi}\leq\phi_{max}\end{split},$$

$$\begin{array}{rl}\underset{\mathbf{\varphi }}{\min }& Z={\tilde{Z}}_1\\ \text{subject to}& \sum _{e\in \Omega }{\rho }^e(\varphi )v^e\le V_{max}\\ & \mathbf{R}(\mathbf{\varphi },\mathbf{u})=0\text{with}\bar{\mathbf{f}}={\bar{\mathbf{f}}}_1+{\bar{\mathbf{f}}}_2\\ & {\varphi }_{min}\le \mathbf{\varphi }\le {\varphi }_{max}\end{array},$$

$$\begin{array}{rl}\underset{\mathbf{\varphi }}{\min }& Z={\tilde{Z}}_2\\ \text{subject to}& \sum _{e\in \Omega }{\rho }^e(\varphi )v^e\le V_{max}\\ & \mathbf{R}(\mathbf{\varphi },{\mathbf{u}}_1)=0\text{with}\bar{\mathbf{f}}={\bar{\mathbf{f}}}_1\\ & \mathbf{R}(\mathbf{\varphi },{\mathbf{u}}_2)=0\text{with}\bar{\mathbf{f}}={\bar{\mathbf{f}}}_2\\ & {\varphi }_{min}\le \mathbf{\varphi }\le {\varphi }_{max}\end{array},$$

$$\begin{array}{r}{\tilde{Z}}_1=\sum _{n=0}^{N_t}({\mathbf{u}}^n{)}^T\mathbf{K}(\varphi ){\mathbf{u}}^n,\end{array}$$

$$\begin{array}{r}{\tilde{Z}}_2=\sum _{n=0}^{N_t}({\mathbf{u}}_1^n{)}^T\mathbf{K}(\varphi ){\mathbf{u}}_1^n+\sum _{n=0}^{N_t}({\mathbf{u}}_2^n{)}^T\mathbf{K}(\varphi ){\mathbf{u}}_2^n,\end{array}$$

$$\begin{split}\tilde{Z}_{1}=\frac{1}{\bar{\xi}}\log_{e}\displaystyle\sum_{n=0}^{N_{t}}\exp\big{(}\bar{\xi}\frac{(\boldsymbol{u}^{n})^{T}\boldsymbol{K}(\phi)\boldsymbol{u}^{n}}{z^{0}_{\mathrm{max}}}\big{)},\end{split}$$

$$\begin{split}\tilde{Z}_{2}=\displaystyle\frac{1}{\bar{\xi}}\log_{e}\sum_{n=0}^{N_{t}}\exp\big{(}\bar{\xi}\frac{(\boldsymbol{u}^{n}_{1})^{T}\boldsymbol{K}(\phi)\boldsymbol{u}^{n}_{1}}{z^{0}_{\mathrm{max}}}\big{)}+\displaystyle\frac{1}{\bar{\xi}}\log_{e}\sum_{n=0}^{N_{t}}\exp\big{(}\bar{\xi}\frac{(\boldsymbol{u}^{n}_{2})^{T}\boldsymbol{K}(\phi)\boldsymbol{u}^{n}_{2}}{z^{0}_{\mathrm{max}}}\big{)},\end{split}$$

$$\begin{array}{r}{\tilde{Z}}_1=({\mathbf{u}}^t{)}^T\mathbf{K}(\varphi ){\mathbf{u}}^t,\end{array}$$

$$\begin{array}{r}{\tilde{Z}}_2=({\mathbf{u}}_1^t{)}^T\mathbf{K}(\varphi ){\mathbf{u}}_1^t+({\mathbf{u}}_2^t{)}^T\mathbf{K}(\varphi ){\mathbf{u}}_2^t,\end{array}$$

Full text

Keywords: Topology optimization, Elastodynamics, Heaviside Projection Method, transient response

Nomenclature

[TABLE]

1 Introduction

Topology optimization as a form-finding methodology for design has developed rapidly during recent years. Its applications cover different fields of engineering such as structural mechanics, fluid flow, heat transfer problems, and so on [4, 5, 27, 2, 21, 1]. For structural problems, it has a demonstrated capability of identifying novel and high performance structures when the operating environment is consistent with assumed load cases used in the optimization. The vast majority of recent studies, however, assume these load cases are static and deterministic. Although significant progress has recently been made considering stochastic dynamic loads described by power spectra, see for example [37], the area of topology optimization under dynamic loads, including the challenging problem of transient loading, is less studied. Some progress has been made in optimization for free vibration problems [9, 7, 31] and forced vibrations [23, 24, 15], including deterministic dynamic and transient loading through time history analysis [16]. This paper presents a gradient-based topology optimization framework for structures subjected to time-dependent and transient loading.

Optimizing the dynamic response of structures has many important applications in different fields of engineering, including structural, mechanical, and aerospace engineering. Analysis is often performed in the time or frequency domains, and topology optimization likewise follows these domains [16, 34]. The objectives of topology optimization in frequency domain are typically minimizing the frequency response of the structure with an excitation frequency or a frequency range (see for example [30]), or minimizing the structural vibration under harmonic excitation (see for example [35, 22]). In the time domain, the objective is often defined as minimizing the strain energy in the structure (maximizing stiffness) within a specific time interval. Due to transient nature of applied loads, designing a structure that preserves the same functionality within different time ranges is a crucial factor in structural dynamics. In the time domain, topology optimization under dynamic loading is fairy difficult and the sensitivities can be costly because many time-dependent analyses are required. Most recent studies, therefore, have focused on simplified methods such as performing topology optimization under a system of equivalent static loads (see for example [6, 17, 20, 14]). The results of topology optimization are used in dynamic analysis of the structure to satisfy the convergence criteria [14]. Model reduction methods using mode displacement and mode acceleration methods have been considered to simplify topology optimization under dynamic loading and reduce computation costs (see for example [36]). Due to the simplifying assumptions of these methods, the engineer may not gain complete control or be able to optimize a specific property of the actual transient response under dynamic loads. For more information about dynamic optimization in time and frequency domains the reader is referred to [33].

Gradient-free optimization methods have been applied to many structural problems with dynamic and transient responses [19, 29, 10]. Unlike gradient-based optimization, these methods do not perform direct sensitivity analysis during the optimization process and instead either generate design alternatives randomly or probabilistically, or approximate design sensitivities using sampling points. These methods are most successful when dimensionality is relatively small, and generally not practical for large scale finite element models with large number of design variables. Gradient-based optimization, on the other hand, offers an effective approach to solve a broad range of topology optimization problems [8]. The computational complexities of sensitivity analysis, however, dramatically increase for transient problems (see for example [25]). Additionally, gradient-based optimization requires an interpolation and penalization scheme to create a continuous representation of material stiffness. The Solid Isotropic Material with Penalization (SIMP) method is the most popular of these methods and uses a power-law interpolation. For more information the reader is referred to [28].

In this paper, we present a topology optimization framework for dynamic response of structures in time domain. The resulting transient topology optimization problem is solved using the gradient-based optimizer Method of Moving Asymptotes (MMA) [32], with time-dependent sensitivities provided via the adjoint method. The SIMP method is used for interpolation of the stiffness of materials. The Heaviside Projection Method (HPM) [12] is used to satisfy the desired minimum length scale of the topological features without adding additional constraints to the problem, thereby influencing manufacturability as well as circumventing well-known instabilities of checkerboard patterns and solution mesh dependency.

2 Physical model

For elastodynamics problems, the time discretized form of the equation of motion is expressed as follows:

[TABLE]

where the superscript $"n"$ denotes current time, $\boldsymbol{M}$, $\boldsymbol{C}$ and $\boldsymbol{K}$ are the mass, viscous damping and stiffness matrices, respectively. In this paper, the viscous damping matrix, $\boldsymbol{C}$, is assumed to be zero. $\boldsymbol{\bar{f}}^{n}$ is the external force vector at current time, $\boldsymbol{\ddot{u}}$, $\boldsymbol{\dot{u}}$ and $\boldsymbol{u}$ are the acceleration, velocity and displacement vectors, respectively. Using the Newmark approach [26], the accelerations and velocities are expressed as follows:

[TABLE]

[TABLE]

where $\tilde{\beta}=0.25$ and $\tilde{\gamma}=0.5$ are parameters for a particular member of the Newmark family [13], and $\Delta t$ is the time step size. Substituting Equations (2) and (3) into (1) yields the following residual equation for the elastodynamics problems:

[TABLE]

This equation is solved for the displacement vector, $\boldsymbol{u}^{n}$, at time $t_{n}$ based on the displacement, $\boldsymbol{u}^{n-1}$, acceleration, $\boldsymbol{\ddot{u}}^{n-1}$, and velocity, $\boldsymbol{\dot{u}}^{n-1}$, vectors at the previous time, $t_{n-1}$. The acceleration, $\boldsymbol{\ddot{u}}^{n}$, and velocity, $\boldsymbol{\dot{u}}^{n}$, are updated recursively using Equations (2) and (3). At the zeroth step, $n=0$, we assume:

[TABLE]

where $\boldsymbol{u}_{\mathrm{int}}$ is the initial condition assumed to be zero.

3 Optimization problem

We present a topology optimization framework for structural problems under transient loading. A representative configuration of an optimization problem considered in this paper is shown in Figure 1. The design domain, $\Omega$, is subjected to a time-dependent loading condition applied at the boundary $\Gamma_{t}$ while displacements are prescribed at the boundary $\Gamma_{u}$. In the optimization examples presented in this paper, we seek to maximize the stiffness (minimize the strain energy) of the structure subjected to transient loading. For transient and dynamic responses, we define the generalized form of the optimization problem in a time-dependent discrete fashion as follows:

@mathmargin

[TABLE]

[TABLE]

where $Z$ is the total objective function, $\boldsymbol{\phi}$ is the vector of design variables, $\boldsymbol{\rho}$ is the vector of densities, $\boldsymbol{u}^{n}$ is the vector of state variables satisfying the discretized governing equations, $\boldsymbol{R}^{n}$, at time steps $n=0,...,N_{t}$. The objective function, $Z$, the equality constraint, $H_{i}$, and the inequality constraint, $G_{j}$, are functions combining performance of the structure at all or select time steps. These functions could simply be a summation to minimize the average transient response, a p-norm function to minimize the maximum response, or some other function to achieve time-dependent goals. The independent design variables are bounded with $\boldsymbol{\phi}_{min}$ and $\boldsymbol{\phi}_{max}$. For compactness, we collect variables at all times to column vectors as follows:

[TABLE]

where $c$ represents either the objective or a constraint. For a transient topology optimization problem, the computational complexities of the sensitivities dramatically increase. The time-dependent adjoint sensitives of a criterion, $\boldsymbol{\tilde{c}}$, the objective or a constraint, with respect to the design variable, $\phi_{k}$, can be written as:

[TABLE]

The total derivative of the state variables with respect to the design variables, $d\boldsymbol{\tilde{u}}/d\phi_{k}$, is computed from the total derivative of the converged residual equation as follows:

[TABLE]

Solving Equation (9) for $d\boldsymbol{\tilde{u}}/d\phi_{k}$ and substituting the result into Equation (8) yields:

[TABLE]

where $\boldsymbol{\tilde{\lambda}}$, vector of adjoint solutions for all time steps, is computed by solving the following adjoint problem:

[TABLE]

For the elastodynamics problem, at the zeroth step, $n=0$, the derivatives of the residual vector with respect to the state vectors, $\partial\boldsymbol{R}^{0}/\partial\boldsymbol{u}^{j}$, are given by:

[TABLE]

where $\boldsymbol{I}$ is the identity matrix. For the $n$th step, $n>0$, the derivatives of $\partial\boldsymbol{R}^{n}/\partial\boldsymbol{u}^{j}$ yield:

[TABLE]

The derivatives of the residual equation with respect to the design variables, $\partial\boldsymbol{\tilde{R}}/\partial\phi_{k}$, can be computed analytically or numerically. From the residual equation, (4), the analytical derivative of the residual with respect to the design variables is described as follows:

[TABLE]

For problems considered here, $\partial\boldsymbol{C}/\partial\phi_{k}=\boldsymbol{0}$. Using the SIMP approach [3], the stiffness of the structure is related to the topology and the Young’s modulus of an element is described as follows:

[TABLE]

where $\rho^{e}$ is the elemental density, $\eta\geq 1$ is the exponent penalty term, $E_{0}^{e}$ is the Young’s modulus of a pure solid material and $\rho_{min}^{e}$ is a small positive number to maintain positive definiteness of the global stiffness matrix. In the HPM, the design variables associated with a material phase are projected onto the finite elements by a regularized Heaviside function [12]. The connection between the design variable space, where the optimization is performed, and the finite element space, where the physical equilibrium is solved, is established through a radial projection. In this method, the projection radius can easily be chosen as the minimum desired feature size, $r_{min}$. The design variables are mapped onto the elements by computing the weighted average, $\mu^{e}(\phi)$, using linear filtering scheme as follows:

[TABLE]

where $x_{j}$ is the position of node $j$, and $x^{e}$ is the central position of element $e$. The elemental density is related to the design variables, $\phi$, as follows:

[TABLE]

where $\beta\geq 0$ dictates the curvature of the regularization. For more information about the HPM the reader is referred to [12, 11].

4 Numerical examples

The method is applied to a structural problem using variety of design considerations. To get insight into the performance of optimization framework, we consider both static and dynamic optimization problems. Throughout this section, a 2D cantilever beam subjected to static and dynamic loading conditions is considered. The schematic of the beam with applied load and boundary conditions is shown in Figure 2, where the loads are applied at points $A$ and $B$. The model and material parameters used in the static and dynamic designs are given in Table 1. We study the influence of static and dynamic loading on the optimized design through numerical examples. Two different optimization cases are considered in static and dynamic designs. In the first case, the optimization problem is solved to minimize the strain energy of the structure subjected to simultaneously applied forces $\boldsymbol{\bar{f}}_{1}=f_{\mathrm{max}}\sin(t)$ and $\boldsymbol{\bar{f}}_{2}=f_{\mathrm{max}}\cos(t)$. While the second optimization case considers minimizing the strain energy in the structure where $\boldsymbol{\bar{f}}_{1}$ and $\boldsymbol{\bar{f}}_{2}$ are applied as independent load cases. More specifically, this case considers two independent forward and sensitivity analyses associated with each $\boldsymbol{\bar{f}}_{1}$ and $\boldsymbol{\bar{f}}_{2}$. For dynamic examples, we further investigate the evolution of the design subjected to dynamic loading by considering three different formulations for the objective function. In these optimization formulations we seek to minimize: (1) sum of the strain energy over the entire transient response, (2) maximum of the strain energy in the time-history response, and (3) the strain energy at a select time. Similar to the static design, we solve two different optimization cases (simultaneously and independently applied forces) for each type of the objective formulation.

The optimization problems are solved by using the MMA [32]. For all optimization examples, we apply a continuation approach on the SIMP exponential penalty term [11]. For numerical modeling, the design domain is discretized with 270 $\times$ 90 elements, and the physical response is predicted with the plane stress condition. The system of structural dynamic equations is solved using the modified Newton-Raphson iteration with an implicit time integration scheme developed by [13].

4.1 Design subjected to static loading

In this section, we explore the characteristics of the proposed optimization framework through a design subjected to static loading. To this end, the 2D cantilever beam, shown in Figure 2, subjected to two different loading conditions is considered. For the static design, the following optimization problems are defined. In the first optimization problem, we seek to:

[TABLE]

where $\boldsymbol{u}$ is the solution corresponding to the structural response subjected to $f_{\mathrm{max}}$ applied at point $A$ and $B$, $v^{e}$ is the elemental volume, and $V_{max}$ is the maximum allowable volume fraction. In the second case, the objective is defined as sum of the strain energy corresponding to sum of the responses of two structures subjected to $f_{\mathrm{max}}$ at point $A$ and $B$, respectively:

[TABLE]

where $\boldsymbol{u}_{1}$ and $\boldsymbol{u}_{2}$ are the solutions corresponding to structural responses subjected to $f_{\mathrm{max}}$ at point $A$ and $B$, respectively. The optimized structures for the static designs are shown in Figure 3. As expected, the optimized structures are different for two optimization problems. Interestingly, in the first case where the structure is subjected to simultaneously applied loads, the shear stresses are zero from the root to midspan (point $B$) leading to the absence of internal bracing, whereas for case (2) the shear stresses are non-zero in this region for both independently applied loads.

4.2 Design subjected to dynamic loading

For the structure subjected to dynamic loading, we consider three different type of objective functions. Similar to the static design, two different cases are considered for each objective function. In the first case, we solve the optimization problem for the structure subjected to simultaneously applied force $\boldsymbol{\bar{f}}_{1}$ and $\boldsymbol{\bar{f}}_{2}$. While in the second case, the optimization problem is solved for the structure subjected to two independent forces $\boldsymbol{\bar{f}}_{1}$ and $\boldsymbol{\bar{f}}_{2}$. The optimization problems are formulated as follows:

case (1):

[TABLE]

case (2):

[TABLE]

where $\tilde{Z}_{1}$ and $\tilde{Z}_{2}$ represent the objective functions. The characteristics of designs subjected to dynamic loading with different objective functions are explored here.

4.2.1 Minimizing the total strain energy

For this optimization problem, we seek to minimize the sum of strain energy over the entire transient response of the structure. Similar to the static design, $\tilde{Z}_{1}$ and $\tilde{Z}_{2}$ are defined as:

[TABLE]

[TABLE]

where $\boldsymbol{u}^{n}$, $\boldsymbol{u}^{n}_{1}$, and $\boldsymbol{u}^{n}_{2}$ are the solution vectors at time step $n$ corresponding to $\boldsymbol{\bar{f}}_{1}+\boldsymbol{\bar{f}}_{2}$, $\boldsymbol{\bar{f}}_{1}$, and $\boldsymbol{\bar{f}}_{2}$, respectively. The optimized structures are shown in Figure 4. The results show two different designs for case (1) and (2). Moreover, both designs differ from the static designs. This is related to the dynamic response of the structure and type of the objective function that considers the structural response at all time steps. Comparison of the optimized strain energy profiles for two cases shows that in case (2), the strain energy associated with the response of structure under $\boldsymbol{\bar{f}}_{2}$ is minor and the design evolves to minimize the governing strain energy.

4.2.2 Minimizing the maximum of strain energy

For this case, we wish to minimize the maximum value of the strain energy over the entire response history of the structure. To construct a differentiable objective function for a min-max optimization problem, we employ the Kreisselmeier-Steinhauser (KS) formulation [18] and formulate the objective functions in the following settings:

[TABLE]

[TABLE]

where $\bar{\xi}$ is a scalar parameter in the KS function which caries between 5 and 200  [18] and $z^{0}_{\mathrm{max}}$ is the maximum value of the strain energy in the time history of the initial design. In this study we set $\bar{\xi}=20$. The optimized designs are shown in Figure 5. The trend of designs is similar to the designs presented in Section 4.2.1; however, the optimized structures are different in case (1) and (2). Moreover, similar to case (2) in Section 4.2.1, the strain energy associated with the response of structure under $\boldsymbol{\bar{f}}_{2}$ is negligible during the optimization. These observations indicate the importance in the timing, location, and magnitude of the applied load for the structural dynamics.

4.2.3 Minimizing strain energy at a designer selected time

Finally, in this section we seek to minimize the strain energy at a time in the time history profile that is selected by the designer. To this end, we define the following objective functions:

[TABLE]

[TABLE]

where $t$ is the designer selected time (e.g. $t=2\pi/5s$, 9th time step) and $\boldsymbol{u}^{t}$ is the solution vector at the selected time. The optimized structures for this type of objective function are given in Figure 6. Unlike previous designs, the optimized designs evolve in a way that minimize the strain energy at the given time. The comparison of the strain energy profiles shows that the optimizer minimizes the strain energy at the given time, however, due to the nature of the defined objective function the strain energy after time $t$ is ignored by the optimizer. Of course this may be impractical, but here we simply wish to demonstrate the ability to influence the behavior of structure at a select time. This is shown in the strain energy profiles given in Figure 6. Moreover, the comparison of the performance of optimized designs under dynamic loading, via the strain energy profiles, shows the importance in the formulation of the objective function and timing, position, and magnitude of applied loads that need to be considered in the structural dynamics analysis, see Figure 7.

5 Conclusions

The gradient-based topology optimization framework was presented to minimize the strain energy in structures subjected to dynamic loading. The time-dependent sensitivities were provided via the adjoint method. The SIMP method was used to interpolate the stiffness of materials and HPM was used to improve the manufacturability and satisfy the minimum length scale in topology-optimized designs. Both static and dynamic optimization problems were considered in this paper. The optimized designs showed significant differences between structures designed for static and dynamic applied loads. Moreover, the optimized dynamic designs revealed the importance in the formulation of the objective function, timing and position of applied loads need to be considered in the structural dynamics design process.

Acknowledgments

This work was supported in part by the National Science Foundation under Grant No. 1538367. This support is gratefully acknowledged.