**1. Introduction**

The demand for creative and optimized designs that cannot be realized using conventional manufacturing techniques has been increasing in recent years owing to the development of additive manufacturing. The design of additive manufacturing (DfAM) has received widespread attention due to developments in the industrial field [1,2]. DfAM is an advanced technique compared to design for manufacturing (DfM) and can be used to overcome the limitations faced by the DfM process. DfAM can be used to produce more creative and high-quality designs since additive manufacturing can produce highly rigid and lightweight products with complicated internal structures. Among the various DfAM techniques, the most widely analyzed is topology optimization, which presents an excellent synergistic effect when combined with additive manufacturing [1–5].

Topology optimization is a scheme which determines the means of the distribution of materials in order to find an optimum topological layout based on the designer's purpose using various simulations [6–8]. The basic concept of topology optimization is to distinguish between the necessary and unnecessary parts in order to best meet the optimization formulation within a given design domain. It has been used in various industries and in research owing to its theoretical simplicity and ability to produce a concise, esthetic, and creative design [9–13]. Despite the several studies conducted on topology optimization, most have been focused on a single part [6–13]. However, in the practical engineering field, each part interacts with its adjacent parts, and most studies have replaced the interaction effects by imposing an equivalent interface condition to address this issue [14–17]. Unfortunately, it is incompatible and impossible in highly nonlinear systems such as contact and impact conditions.

In order to consider the effects of adjacent parts, Strömberg et al. suggested a topology optimization framework which simultaneously considers multiple deformable parts

through a node-to-node contact condition [18]. However, this is only feasible in the matching interface meshes owing to the features of the node-to-node contact condition. Jeong et al. presented a new topology optimization technique based on a surface-to-surface contact condition to handle both the matching and non-matching interface meshes [19]. The surface-to-surface contact condition was treated by the condensed mortar method and a filtering scheme considering the size of the element was introduced to alleviate some discontinuities present on the non-matching interface. Recently, Fernandez et al. proposed a topology optimization technique to handle multiple deformable bodies in static and large deformation contacts based on a mortar formulation approach [20].

This paper presents a method to determine an optimum topological layout for the linear dynamic impact and frictionless contact conditions to present a more generalized application of topology optimization. Firstly, the condensed mortar method is used to handle the non-matching interface, which is inevitably caused by the relative motion near the impact region. The condensed mortar method presents several advantages when applied in topology optimization because it can be applied in the same manner as that of conventional structural analysis, even in the presence of a non-matching interface [19]. Since previous studies only considered a linear static problem [19], the proposed approach presents a way of applying the condensed mortar method in a linear dynamic topology optimization. Secondly, the solid isotropic method with penalization (SIMP) method is used to determine the optimum topological layout [13]. The mean dynamic compliance and mean squared displacement of a target DOF are considered as an objective function in this study to obtain an optimum topological layout for a linear dynamic system [13]. Furthermore, a filtering scheme that considers the size of an element is applied to handle discontinuities along the non-matching interface [19]. As a result, the approach proposed in this paper can be easily used in various applications including analysis and design because the non-matching interface are simply handled by applying the condensed localized mortar method to the well-known linear dynamic analysis methods.

The remainder of this paper is organized as follows. Section 2 briefly introduces the concept of the condensed mortar method, which is a scheme used to handle non-matching interfaces. The validity of the condensed mortar method is verified through a simple verification example for the linear dynamic analysis. Section 3 describes a method to determine an optimum topological layout in the linear dynamic problems. The condensed mortar method is applied to handle the non-matching interface on the impact regions, based on the conventional dynamic response topology optimization framework. Section 4 presents the verification of the validity of the proposed method using an impact problem in which a beam and a flying block are impacted by frictionless contact conditions. Lastly, Section 5 presents the concluding remarks.

#### **2. Linear Dynamic Analysis with Condensed Localized Mortar Method**

This section presents a linear dynamic analysis framework using the condensed localized mortar method to handle the impact and frictionless contact conditions. Figure 1 illustrates the overall analysis process.

When the analysis is initiated, it is necessary to define the system and to determine how the system is considered at the problem start step. Subsequently, the contact region between the adjoined parts is searched in the case of the impact and frictionless contact problem. A linear dynamic structural analysis considering the contact conditions is performed in the event of an impact or contact; otherwise, a linear dynamic structural analysis is performed for one part. In this paper, only the impact and frictionless contact conditions are considered among the various contact conditions (i.e., interface-tie, friction or frictionless contact, and impact conditions). Additionally, the contact conditions are handled using the condensed mortar method. A detailed explanation of this is presented in the following subsections.

The overall process ends when the time *t* is greater than the final time *tf* . Otherwise, the variables including the time and information for node, displacement, velocity, and acceleration are updated, and the process is repeated, as shown in Figure 1.

**Figure 1.** Flowchart for linear dynamic contact analysis.

#### *2.1. Linear Dynamic Finite Element Analysis Considering Contact Condition*

This subsection presents the basic concept of the 2D linear dynamic finite element analysis used to solve the impact and frictionless contact problems by employing the condensed mortar method. Figure 2 presents a schematic diagram of the 2D contact problem. Even though only two contacting parts are considered in Figure 2, it can be easily extended to multiple parts. The initial boundary value problem of the general structural mechanics comprises a set of combined second-order partial differential equations and a given set of initial and boundary conditions [21].

**Figure 2.** Schematic description for contact problem.

The initial boundary value problem for the linear dynamic structural mechanics is given as:

$$\begin{array}{c} \text{Div}\mathbf{P} + \mathbf{b} = \rho \ddot{\mathbf{u}} \text{ in } \Omega \times [0, \text{ T}]\\ \mathbf{u} = \stackrel{\leftarrow}{\mathbf{u}} \text{ on } \Gamma\_{\boldsymbol{u}} \times [0, \text{ T}]\\ \mathbf{P} \cdot \mathbf{N} = \stackrel{\leftarrow}{\mathbf{t}} \text{ on } \Gamma\_{\boldsymbol{t}} \times [0, \text{ T}] \end{array} \tag{1}$$

where Ω is the domain occupied by all the material points, Γ*<sup>u</sup>* and Γ*<sup>t</sup>* denote the Dirichlet

and Neumann boundary, where displacements ¯ **<sup>u</sup>** and tractions ¯ **t** are prescribed. **P** represents the first Piola-Kirchhoff stress tensor, **b** represents an external body force vector, *ρ* is the density, **u** is the displacement vector, and the overdot indicates the time derivative.

In order to address the impact condition with the frictionless contact between adjoined parts, the Karush-Kuhn-Tucker conditions must be satisfied in the normal direction

$$\emptyset \mathfrak{g}\_n \ge 0, \ \mathfrak{r}\_n \le 0, \ \mathfrak{r}\_n \mathfrak{g}\_n = 0 \tag{2}$$

where *τ<sup>n</sup>* represents the normal component of the contact traction and *gn* indicates a gap function in the normal direction, which is defined by:

$$\mathbf{g}\_n = -\mathbf{n} \,\mathbf{x}^{(1)} \left[ \mathbf{x}^{(1)} - \tilde{\mathbf{x}}^{(2)} \right] \tag{3}$$

where **<sup>~</sup> x** (2) is defined by projecting the point **x**(1) onto the surface Γ*<sup>c</sup>* in the domain Ω<sup>2</sup> along the normal direction **<sup>n</sup>**. The projected point **<sup>~</sup> x** (2) can be obtained by using a well-known contact search algorithm (e.g., [22]).

When the finite element method is employed to perform structural analysis based on Equations (1)–(3), the momentum equation under an external load can be given by:

$$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{f}\_{\text{int}}(\mathbf{u}) + \mathbf{f}\_{\varepsilon}(\mathbf{u}, \lambda) - \mathbf{f}\_{\text{ext}} = 0 \tag{4}$$

where **M** is the global mass matrix and **f**int, **f**ext and **f***<sup>c</sup>* represent the global internal, external, and contact force vectors, respectively. Here, the mass matrix **M** is assumed to be constant.

In order to obtain the discrete single unknown **u***n*+<sup>1</sup> in single step '*n* + 1 , Newmark approximation [23] is used to describe other quantities as a function of **u***n*+1.

$$\begin{array}{l} \dot{\mathbf{u}}\_{n+1} = \frac{\gamma}{\beta \Delta t} (\mathbf{u}\_{n+1} - \mathbf{u}\_{n}) - \frac{\gamma - \beta}{\beta} \dot{\mathbf{u}}\_{n} - \frac{\gamma - 2\beta}{2\beta} \Delta t \ddot{\mathbf{u}}\_{n} \\\ \ddot{\mathbf{u}}\_{n+1} = \frac{\gamma}{\beta \Delta t^{2}} (\mathbf{u}\_{n+1} - \mathbf{u}\_{n}) - \frac{1}{\beta \Delta t} \dot{\mathbf{u}}\_{n} - \frac{1 - 2\beta}{2\beta} \ddot{\mathbf{u}}\_{n} \end{array} \tag{5}$$

where *β* and *γ* represent the algorithmic parameters for the Newmark method, typically 2*β* = *γ* = 1/2.

#### *2.2. Application of Condensed Localized Mortar Method in Linear Dynamic Analysis*

Based on the basic concept and framework of the linear dynamic analysis considering frictionless contact conditions, this subsection presents a method to apply the condensed localized mortar method [19] in the linear dynamic analysis to derive the contact force vector **f***c* in Equation (4).

The mortar method is the most widely used approach to solve the interface, contact, impact, and multi-physics problems by introducing Lagrange multipliers to impose the interface compatibility condition in a weak sense [24–27]. The mortar methods can be divided into classical and localized versions depending on how the interface compatibility conditions and Lagrange multipliers are defined. The classical mortar method is the most widely used approach for treating the interface conditions by introducing Lagrange multipliers as the interface pressure, as shown in Figure 3b. A fictitious frame is introduced in the localized mortar method, and the interface compatibility condition is imposed through the frame domain, as shown in Figure 3c [19,28–32].

**Figure 3.** Simple patch test problem: (**a**) a continuum structure composed of two subdomains, and each subdomain is connected using the (**b**) classical and (**c**) localized version.

This study employs the localized mortar method to impose the interface compatibility conditions because it is successfully utilized in various problems by uniquely defining the interface compatibility conditions through the frame domain [19,28–32]. In the case of the localized mortar method, the interface compatibility condition between the newly introduced frame domain and the adjacent two parts, as shown in Figure 3c, is given by

$$\mathbf{u}^{(1)} - \mathbf{u}^{(2)} = 0 \implies \begin{cases} \mathbf{u}^{(1)} - \mathbf{u}^{(f)} = 0 \\ \mathbf{u}^{(2)} - \mathbf{u}^{(f)} = 0 \end{cases} \tag{6}$$

Subsequently, two sets of Lagrange multipliers, i.e., λ(1) and λ(2) , are defined along the interface to impose the interface capability condition (6) in a weak sense. By substituting the interface condition (6) and two sets of Lagrange multipliers into the momentum equation of the mortar method (4) and organizing the matrix form, the total system equation of the localized mortar method in a linear dynamic condition is given by:

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ **M**(1) 0 000 0 **M**(2) 000 0 0 000 0 0 000 0 0 000 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ **.. <sup>u</sup>**(1) **.. <sup>u</sup>**(2) **..** λ (1) **..** λ (2) **.. <sup>u</sup>**(*f*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ **K**(1) 0 **B**1**Q**<sup>1</sup> 0 0 0 **K**(2) 0 **B**2**Q**<sup>2</sup> 0 **Q***<sup>T</sup>* <sup>1</sup> **B***<sup>T</sup>* <sup>1</sup> 0 00 −**W**<sup>1</sup> 0 **Q***<sup>T</sup>* <sup>2</sup> **B***<sup>T</sup>* <sup>2</sup> 0 0 −**W**<sup>2</sup> 0 0 <sup>−</sup>**W***<sup>T</sup>* <sup>1</sup> <sup>−</sup>**W***<sup>T</sup>* <sup>2</sup> 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ **u**(1) **u**(2) λ(1) λ(2) **u**(*f*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ **f** (1) **f** (2) **G**(1) **G**(2) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , **Q**<sup>1</sup> = **N***<sup>T</sup> u*Γ<sup>1</sup> **N***λ*1*d*Γ, **Q**<sup>2</sup> = **N***<sup>T</sup> <sup>u</sup>*Γ2**N***λ*2*d*Γ, **W**<sup>1</sup> = **N***<sup>T</sup> <sup>λ</sup>*1**N***u f d*Γ, **W**<sup>2</sup> = **N***<sup>T</sup> <sup>λ</sup>*2**N***u f d*Γ (7)

Γ*I* Γ*I* Γ*I* Γ*I* **G**(1) = Γ*I* **N***<sup>T</sup> <sup>λ</sup>*1*g*<sup>1</sup> *<sup>f</sup> <sup>d</sup>*Γ, **<sup>G</sup>**(2) <sup>=</sup> Γ*I* **N***<sup>T</sup> <sup>λ</sup>*2*g*<sup>2</sup> *<sup>f</sup> d*Γ

> where **B**<sup>1</sup> and **B**<sup>2</sup> are the Boolean matrices representing the boundary displacement components, and **Q**1, **Q**2, **W**1, and **W**<sup>2</sup> are the projection matrices, which represent the interface condition between the interface displacements and frame displacement. Additionally, **G**(1) and **G**(2) indicate the non-penetration condition caused by the KKT condition (2) and the interface compatibility condition (6).

> The condensed localized mortar method, which was developed to be easily applied in topology optimization algorithms even in the presence of a non-matching interface, is employed [19]. This can be obtained by condensing the Lagrange multipliers and overlapping interface displacements using the localized mortar method (7).

> When the total system equation is condensed, the numerical efficiency is increased drastically, and the total system equation can be used universally because it can be considered equivalent to the framework of the general structural finite element analysis [19].

> Before proceeding with the condensed localized mortar method, the displacement **u**(*n*) at each part must be divided into the internal component **<sup>u</sup>**(*n*) *<sup>i</sup>* and the interface component

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

**<sup>u</sup>**(*n*) *<sup>b</sup>* as shown in Figure 4. When the displacements are decomposed, the total system Equation (7) for the localized mortar method is rewritten as:

**Figure 4.** Division of the total displacements into the internal and interface components.

The following two steps must be performed to derive the condensed localized mortar method. Firstly, two relations must be obtained through the fifth and sixth rows of the total system Equation (8), to condense the boundary displacements **<sup>u</sup>**(1) *<sup>b</sup>* , **<sup>u</sup>**(2) *<sup>b</sup>* into the frame displacement **u**(*f*).

$$\begin{aligned} \text{5th- row}: \mathbf{Q}^{(1)T} \mathbf{u}\_b^{(1)} - \mathbf{W}^{(1)} \mathbf{u}^{(f)} &= 0 \rightarrow \mathbf{u}\_b^{(1)} = \left(\mathbf{Q}^{(1)T}\right)^{-1} \mathbf{W}^{(1)} \mathbf{u}^{(f)} = \mathbf{C}\_{1f} \mathbf{u}^{(f)} \\ \text{6th- row}: \mathbf{Q}^{(2)T} \mathbf{u}\_b^{(2)} - \mathbf{W}^{(2)} \mathbf{u}^{(f)} &= 0 \rightarrow \mathbf{u}\_b^{(2)} = \left(\mathbf{Q}^{(2)T}\right)^{-1} \mathbf{W}^{(2)} \mathbf{u}^{(f)} = \mathbf{C}\_{2f} \mathbf{u}^{(f)} \end{aligned} \tag{9}$$

Secondly, the second and fourth rows of the total system Equation (8) are used as follows to eliminate the Lagrange multipliers, λ(1) and λ(2) :

$$\begin{array}{c} \text{2nd- row}: \mathbf{M}^{(1)}\_{bi}\ddot{\mathbf{u}}^{(1)}\_{i} + \mathbf{K}^{(1)}\_{bi}\mathbf{u}^{(1)}\_{i} + \mathbf{M}^{(1)}\_{bb}\ddot{\mathbf{u}}^{(1)}\_{b} + \mathbf{K}^{(1)}\_{bb}\mathbf{u}^{(1)}\_{b} + \mathbf{Q}^{(1)}\lambda^{(1)} = \mathbf{f}^{(1)}\_{b} \\ \rightarrow \lambda^{(1)} = \left(\mathbf{Q}^{(1)}\right)^{-1}\left(\mathbf{f}^{(1)}\_{b} - \mathbf{M}^{(1)}\_{bi}\mathbf{u}^{(1)}\_{i} - \mathbf{K}^{(1)}\_{bi}\mathbf{u}^{(1)}\_{i} - \mathbf{M}^{(1)}\_{bi}\mathbf{u}^{(1)}\_{b} - \mathbf{K}^{(1)}\_{bb}\mathbf{C}\_{1f}\mathbf{u}^{(f)}\right) \\ \quad \text{4th- row}: \mathbf{M}^{(2)}\_{bi}\ddot{\mathbf{u}}^{(2)}\_{i} + \mathbf{K}^{(2)}\_{bi}\mathbf{u}^{(2)}\_{i} + \mathbf{M}^{(2)}\_{bb}\ddot{\mathbf{u}}^{(2)}\_{b} + \mathbf{K}^{(2)}\_{bb}\mathbf{u}^{(2)}\_{b} + \mathbf{Q}^{(2)}\lambda^{(2)} = \mathbf{f}^{(2)}\_{b} \\ \rightarrow \lambda^{(2)} = \left(\mathbf{Q}^{(2)}\right)^{-1}\left(\mathbf{f}^{(2)}\_{b} - \mathbf{M}^{(2)}\_{bi}\dot{\mathbf{u}}^{(2)}\_{i} - \mathbf{K}^{(2)}\_{bi}\mathbf{u}^{(2)}\_{i} - \mathbf{M}^{(2)}\_{bb}\dot{\mathbf{u}}^{(2)}\_{b} - \mathbf{K}$$

When the relations (9) and (10) are substituted into the total system Equation (8), the Lagrange multipliers <sup>λ</sup>(1) and <sup>λ</sup>(2) and the overlapped boundary displacements **<sup>u</sup>**(1) *<sup>b</sup>* and **<sup>u</sup>**(2) *<sup>b</sup>* are condensed by the internal displacements **<sup>u</sup>**(1) *<sup>i</sup>* and **<sup>u</sup>**(2) *<sup>i</sup>* and the frame displacement **u**(*f*). Lastly, the matrix form of the total system equation can be obtained as follows by summarizing and re-organizing Equation (8):

**M**ˆ **.. u**ˆ + **K**ˆ **u**ˆ = **F**ˆ where **M**ˆ = ⎡ ⎢ ⎢ ⎣ **<sup>M</sup>**(1) *ii* <sup>0</sup> **<sup>M</sup>**(1) *ib* **C**<sup>1</sup> *<sup>f</sup>* <sup>0</sup> **<sup>M</sup>**(2) *ii* **<sup>M</sup>**(2) *ib* **C**<sup>2</sup> *<sup>f</sup>* **C***T* <sup>1</sup> *<sup>f</sup>* **<sup>M</sup>**(1) *bi* **<sup>C</sup>***<sup>T</sup>* <sup>2</sup> *<sup>f</sup>* **<sup>M</sup>**(2) *bi* **<sup>C</sup>***<sup>T</sup>* <sup>1</sup> *<sup>f</sup>* **<sup>M</sup>**(1) *bb* **<sup>C</sup>**<sup>1</sup> *<sup>f</sup>* <sup>+</sup> **<sup>C</sup>***<sup>T</sup>* <sup>2</sup> *<sup>f</sup>* **<sup>M</sup>**(2) *bb* **C**<sup>2</sup> *<sup>f</sup>* ⎤ ⎥ ⎥ ⎦, **<sup>u</sup>**<sup>ˆ</sup> <sup>=</sup> ⎡ ⎢ ⎣ **<sup>u</sup>**(1) *i* **<sup>u</sup>**(2) *i* **u**(*f*) ⎤ ⎥ ⎦, **.. u**ˆ = ⎡ ⎢ ⎢ ⎣ **.. <sup>u</sup>**(1) *i* **.. <sup>u</sup>**(2) *i* **.. <sup>u</sup>**(*f*) ⎤ ⎥ ⎥ ⎦ **K**ˆ = ⎡ ⎢ ⎢ ⎣ **<sup>K</sup>**(1) *ii* <sup>0</sup> **<sup>K</sup>**(1) *ib* **C**<sup>1</sup> *<sup>f</sup>* <sup>0</sup> **<sup>K</sup>**(2) *ii* **<sup>K</sup>**(2) *ib* **C**<sup>2</sup> *<sup>f</sup>* **C***T* <sup>1</sup> *<sup>f</sup>* **<sup>K</sup>**(1) *bi* **<sup>C</sup>***<sup>T</sup>* <sup>2</sup> *<sup>f</sup>* **<sup>K</sup>**(2) *bi* **<sup>C</sup>***<sup>T</sup>* <sup>1</sup> *<sup>f</sup>* **<sup>K</sup>**(1) *bb* **<sup>C</sup>**<sup>1</sup> *<sup>f</sup>* <sup>+</sup> **<sup>C</sup>***<sup>T</sup>* <sup>2</sup> *<sup>f</sup>* **<sup>K</sup>**(2) *bb* **C**<sup>2</sup> *<sup>f</sup>* ⎤ ⎥ ⎥ ⎦, **<sup>F</sup>**<sup>ˆ</sup> <sup>=</sup> ⎡ ⎢ ⎢ ⎣ **f** (1) *i* **f** (2) *i* **C***T* <sup>1</sup> *<sup>f</sup>***f** (1) *<sup>b</sup>* <sup>+</sup> **<sup>C</sup>***<sup>T</sup>* <sup>2</sup> *<sup>f</sup>***f** (2) *b* ⎤ ⎥ ⎥ ⎦ (11)

The internal displacements **<sup>u</sup>**(1) *<sup>i</sup>* and **<sup>u</sup>**(2) *<sup>i</sup>* and the frame displacement **<sup>u</sup>**(*f*) can be obtained directly by solving the total system equation of the condensed localized mortar method (11). The condensed variables **<sup>u</sup>**(1) *<sup>b</sup>* , **<sup>u</sup>**(2) *<sup>b</sup>* , <sup>λ</sup>(1) and <sup>λ</sup>(2) can then be calculated using the relation (9) and (10), if required.

#### *2.3. Verification Example for Linear Dynamic Impact and Frictionless Contact Problem*

The problem of ring contacts to the rigid support beam [33] is considered to ensure that the condensed localized mortar method works well in the linear dynamic impact and frictionless contact conditions. Figure 5a presents the problem definition.

As shown in Figure 5a, both side of the rigid support with a size of 24 by 3 is completely fixed and the displacement boundary condition is imposed on the top surface of the ring so that the ring goes down constantly. All the parts are modeled with a 2D linear quadrilateral element for both parts. In the verification example, the elastic ring consists of an inner and outer ring with the same thickness and different materials [33]. In order to ensure that the inner ring is stiffer than the outer ring, the material properties of the inner and outer rings are defined as follows:

$$\begin{array}{l}E\_{\text{inner}} = 1E + 05, \nu = 0.3\\E\_{\text{outer}} = 1E + 03, \nu = 0.3\end{array} \tag{12}$$

By imposing the condition of pressing on the top side of the ring, as shown in Figure 5a, the finite deformable ring is in contact with a rigid support, and a non-matching interface is then observed on the interface. The contact condition is solved by using the condensed localized mortar method (11).

A linear dynamic finite element analysis is conducted in 120 steps to solve the verification example. The deformed configurations are presented in Figure 5b–f. The contact area is changed from one spot to one region and then expanded to two regions to distribute the pressure by increasing the contact surface in order to support a given boundary condition as the analysis proceeds. These analysis results demonstrate identical appearances as shown in reference [33].

This result indicates that the condensed localized mortar method used in this study works well even in the linear dynamic impact and frictionless contact conditions, and that the contact regions are changed. Therefore, the interface is treated by the condensed localized mortar method to handle the topology optimization with the impact and frictionless contact conditions.

## **3. Dynamic Response Topology Optimization with Condensed Localized Mortar Method**

Topology optimization is widely applied as a conceptual design that can be used to derive well-matched structures without any intervention from the designer [6–8]. The SIMP method is the most widely used approach among the various topology optimization techniques [9–12]; it is used to gradually find an optimum layout by determining the material existence or absence of each element in a given design domain. The SIMP method uses a relative density function to determine the existence of a material *η*.

$$\eta(\mathbf{x}) = \begin{cases} 1 \text{ if } \mathbf{x} \in \Omega \\ 0 \text{ if } \mathbf{x} \notin \Omega \end{cases} \tag{13}$$

where *x* indicates the location of the finite element in the design domain, Ω. From Equation (13), it can be inferred that the material exists when the relative density function *η* is 1, and the material does not exist when the relative density function *η*, is 0. In the modified SIMP method [13], material properties such as Young's modulus and density of each element are defined based on the function of relative density *η* as given below:

$$\begin{aligned} E\_{\mathfrak{c}} &= E\_0 \left( a \eta\_{\mathfrak{c}}^p + (1 - a) \eta\_{\mathfrak{c}} \right) \\ \rho\_{\mathfrak{c}} &= \rho\_0 \eta\_{\mathfrak{c}} \end{aligned} \tag{14}$$

where *E*<sup>0</sup> and *ρ*<sup>0</sup> are the Young's modulus and density of the material, respectively, *α* is a positive constant that is typically defined as 15/16, and *p* is a penalization parameter that is typically set as 3.

#### *3.1. Optimization Formulation with Condensed Localized Mortar Method*

The optimization formulation used to obtain an optimum topological layout in a linear dynamic system [13,34,35] is represented by Equation (15), where the condensed localized mortar method described in the previous section is adopted by the matrices **M**ˆ , **K**ˆ , and vectors **u**ˆ(*t*), **.. u**ˆ(*t*) and **F**ˆ(*t*) to handle the impact problem with the non-matching interface.

$$\begin{aligned} \underset{\eta}{\text{minimize}} \quad & f\_d \\ \text{subject to} \begin{cases} \quad & \mathbf{\tilde{M}} \ddot{\mathbf{u}}(t) + \mathbf{\tilde{K}u}(t) = \mathbf{f}(t) \; \forall t \in [0, t\_f] \\ \quad & g(\eta) = V(\eta) - V\_{\text{re}\eta} = \sum\_{c=1}^{N} \eta\_c v\_c - V\_{\text{re}\eta} \le 0 \\ \quad & 0 < \eta\_{\text{min}} \le \eta \le 1 \end{cases} \end{aligned} \tag{15}$$

The optimization formulation (15) is organized to minimize the objective function *Jd* using the momentum equation and volume constraint, *g*(*n*) [13]. In this study, the objective function *Jd* is defined for mean dynamic compliance to find a stiff structure as follows:

$$J\_d = \frac{1}{t\_f} \int\_0^{t\_f} \mathfrak{c}(\hat{\mathfrak{u}}(t), \mathfrak{n}) \mathrm{d}t, \mathfrak{c}(\hat{\mathfrak{u}}(t), \mathfrak{n}) = \mathfrak{F}^T(t)\hat{\mathfrak{u}}(t) \text{ where } \frac{\partial \mathfrak{c}}{\partial \eta\_{\mathfrak{c}}} = 0, \frac{\partial \mathfrak{c}}{\partial \hat{\mathfrak{u}}} = \mathfrak{F}(t) \tag{16}$$

Additionally, for the mean squared displacement of a target, DOF is given as:

$$J\_d = \frac{1}{t\_f} \int\_0^{t\_f} c(\hat{\mathfrak{u}}(t), \mathfrak{n}) \mathrm{d}t, c(\hat{\mathfrak{u}}(t), \mathfrak{n}) = (\mathrm{L}\hat{\mathfrak{u}}(t))^2 \text{ where } \frac{\partial c}{\partial \eta\_\varepsilon} = 0, \frac{\partial c}{\partial \hat{\mathfrak{u}}} = 2\mathfrak{u}(t)\mathrm{L} \tag{17}$$

where **L** represents the role in selecting the target DOF from the displacement **u**ˆ(*t*).

The optimum topological layout that minimizes the objective function can be determined, which includes the mean dynamic compliance and mean squared displacement of a target degree-of-freedom while satisfying the constraints, by performing the optimization formulation defined in Equations (15)–(17).

The sensitivity of the optimization formulation (15) must be calculated when a gradient-based optimization algorithm is used. Since the topology optimization presents several design variables by defining the relative density function *η*, for each element, the adjoint method is more appropriate than the direct method [6,7,13]. The adjoint variables λadj must be introduced to simplify the objective function in order to use the adjoint method to calculate the sensitivities of the objective function. When the adjoint variable λadj is inserted with the additional adjoint equation **<sup>M</sup>**<sup>ˆ</sup> **.. <sup>u</sup>**<sup>ˆ</sup> <sup>+</sup> **<sup>K</sup>**<sup>ˆ</sup> **<sup>u</sup>**<sup>ˆ</sup> <sup>−</sup> **<sup>F</sup>**<sup>ˆ</sup> the objective function *Jd* can be rewritten as:

$$J\_d \cong J\_d + \lambda^T \int\_0^{t\_f} (\mathbf{\hat{M}}\ddot{\mathbf{u}} + \mathbf{\hat{K}}\dot{\mathbf{u}} - \mathbf{f}) \, \text{d}t \tag{18}$$

Since the additional adjoint equation **<sup>M</sup>**<sup>ˆ</sup> **.. <sup>u</sup>**<sup>ˆ</sup> <sup>+</sup> **<sup>K</sup>**<sup>ˆ</sup> **<sup>u</sup>**<sup>ˆ</sup> <sup>−</sup> **<sup>F</sup>**<sup>ˆ</sup> is always zero, as illustrated in Equation (11), it can be confirmed that the objective function *Jd* does not change even when the adjoint variable and the adjoint equation are added.

The sensitivity of the objective function is calculated by using the adjoint method, as follows:

$$\frac{\partial l\_d}{\partial \eta\_\varepsilon} = \frac{\partial l\_d}{\partial \eta\_\varepsilon} + \int\_0^{t\_f} \lambda^T \frac{\partial}{\partial \eta\_\varepsilon} (\mathbf{M}\tilde{\mathbf{u}} + \mathbf{K}\hat{\mathbf{u}} - \mathbf{F}) \mathrm{d}t = \int\_0^{t\_f} \left[ \left\{ \frac{\partial c}{\partial \eta\_\varepsilon} + \frac{\partial \tilde{\mathbf{u}}}{\partial \eta\_\varepsilon} \, \frac{\partial \mathbf{c}}{\partial \tilde{\mathbf{u}}} \right\} + \left\{ \lambda^T \frac{\partial}{\partial \eta\_\varepsilon} (\mathbf{M}\tilde{\mathbf{u}} + \mathbf{K}\hat{\mathbf{u}}) \right\} \right] \mathrm{d}t \tag{19}$$

Subsequently, by integrating-by-parts and re-organizing, Equation (19) yields:

$$\frac{\partial l\_{d}}{\partial \boldsymbol{\eta}\_{\varepsilon}} = \int\_{0}^{t\_{f}} \left[ \boldsymbol{\lambda}^{T} \frac{\partial \mathbf{\hat{M}}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \boldsymbol{\tilde{\mathbf{u}}} + \boldsymbol{\lambda}^{T} \frac{\partial \mathbf{\hat{K}}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \boldsymbol{\hat{u}} + \frac{\partial \boldsymbol{c}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \right] d\mathbf{t} + \int\_{0}^{t\_{f}} \left[ \frac{\partial \mathbf{\hat{u}}^{T}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \left( \frac{\partial \mathbf{c}}{\partial \boldsymbol{\tilde{u}}} + \boldsymbol{\hat{M}} \boldsymbol{\tilde{\lambda}} + \boldsymbol{\hat{K}} \boldsymbol{\lambda} \right) \right] d\mathbf{t} + \left. \left[ \frac{\partial \mathbf{\hat{u}}^{T}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \boldsymbol{\hat{M}} \lambda - \frac{\partial \mathbf{\hat{u}}^{T}}{\partial \boldsymbol{\eta}\_{\varepsilon}} \boldsymbol{\hat{M}} \boldsymbol{\hat{\lambda}} \right] \right|\_{t=t\_{f}} \tag{20}$$

In Equation (20), to obtain the adjoint variable λadj that is arbitrarily inserted to simplify the sensitivity process, the following differential equation must be solved:

$$\mathbf{\hat{M}}\ddot{\mathbf{\lambda}} + \mathbf{\hat{K}}\lambda = -\frac{\partial \mathbf{c}}{\partial \mathbf{\hat{u}}'} \ t \in \left[0, t\_f\right] \text{ where } \lambda\_{t\_f} = 0, \dot{\lambda}\_{t\_f} = 0 \tag{21}$$

where Equation (21) is a similar form of the equation of motion (11), which is widely known. Therefore, the adjoint variable is transformed by Λ(*s*) = λadj(*tf* − *s*) in order to make an equivalent form, and Equation (21) is re-written as:

$$\mathbf{\hat{M}}\dot{\boldsymbol{\Lambda}}(s) + \mathbf{\hat{K}}\boldsymbol{\Lambda}(s) = -\frac{\partial \boldsymbol{c}}{\partial \mathbf{\hat{u}}}\bigg|\_{t\_f - s}, \text{ s} \in \left[0, t\_f\right] \text{ where } \boldsymbol{\Lambda}(0) = 0, \dot{\boldsymbol{\Lambda}}(0) = 0 \tag{22}$$

Since the form of Equation (22) is exactly identical to that of Equation (11), the transformed adjoint variable Λ(*s*) can be easily calculated, along with the adjoint variable λadj through an inverse transformation.

#### *3.2. Verification Example for Dynamic Response Topology Optimization*

A Michell-type structure problem is used to validate the optimization formulation (15) for the linear dynamic system. The geometry and boundary conditions are shown in Figure 6a,b. The material properties of the entire design domain are assumed to be those of structural steel (*<sup>E</sup>* = <sup>2</sup> × 1011, Poisson's ratio is 0.3, and density *<sup>ρ</sup>* = 7890). The total time step is considered as 50.

**Figure 6.** Problem definition for Michell-type structure problem: (**a**) schematic description (**b**) load condition.

Topology optimization is performed to obtain the stiffest design while satisfying 50% of the volume using the optimization formulation (15) and (16). The adjoint method is used for sensitivity analysis. In order to identify the differences and effectiveness of implementation in the linear dynamic system before performing the dynamic response topology optimization, the static response topology optimization under the same conditions and optimization formulations is performed as shown in Figure 7.

In order to perform topology optimization, it is necessary to confirm whether the dynamic analysis is performed adequately due to which, the Michell-type structure is analyzed first. Figure 8 shows the vertical directional displacements at point "A". The red color represents the result of the static analysis, and the other results indicate the displacement of the dynamic analysis corresponding to the final time *tf* . In the verification example, analysis and optimization is conducted with the final time *tf* set to 0.01, 0.05, 0.1, and 0.5 to check the response in the topology optimization for the dynamic system.

**Figure 7.** Optimized results for Michell-type structure with static analysis: (**a**) Optimized layout (**b**) convergence history.

F'LVSODFHPHQW DW³\$´ ZKHQݐ = 0.1 G'LVSODFHPHQW DW³\$´ZKHQݐ = 0.

**Figure 8.** Displacement histories of the tip point for Michell-type structure.

As shown in Figure 8, the results of the displacement at point "A" are converged to the quasi-static results as the final time *tf* increases. Figure 8a shows the displacement at point "A" when the final time *tf* is 0.01 s, and the structure acts as an impact state because the load is applied for a very short period of time. This indicates that the dynamic effect produced by the impact phenomenon is essential in the state of Figure 8a with a final time *tf* of 0.01 s.

In the case of the impact state as shown in Figure 9a, the material is concentrated near the point "A" where the load is applied due to the effect of inertia. Conversely, in the case of quasi-static cases, as shown in Figure 9b–d, the optimized topological layout is identical to the result of the static analysis presented in Figure 7 because the structure acts as a quasi-static state when the final time *tf* is greater than 0.05.

The verification example is used to verify that the dynamic response topology optimization adequately derives the optimized layout when the dynamic effect is critical. As shown in Figure 9, the material is concentrated near the impact region due to the inertia effect in the impact state. However, in a quasi-static state, the optimized topological layout is obtained as the result of the static analysis as the final time *tf* increases. This indicates that the dynamic response topology optimization is applied well in both the impact and quasi-static states when the formulation (15) is used.

## **4. Numerical Example: Impact between Fixed Beam and Flying Block**

Hitherto, two methods have been employed to find an optimum topological layout in the linear dynamic impact problem: the condensed localized mortar method and the dynamic response topology optimization approach. The condensed localized mortar method used to handle the non-matching interface at the impact and frictionless contact region was described in Section 2, and the simple frictionless contact problem is solved to verify its effectiveness in the impact state. The dynamic response topology optimization was introduced in Section 3, and an optimization formulation with the condensed mortar method was presented. Additionally, by performing topology optimization on the verification example, the effectiveness of the dynamic response topology optimization is confirmed when the beam is in the impact state.

In this section, topology optimization is performed by assuming a situation in which two independently behaving objects are in the impact and frictionless contact conditions to examine the applicability of the proposed approach to the linear dynamic impact problem. The problem definition is presented in Figure 10. A block is separated at 0.05 intervals on a lower beam which is completely fixed on both sides, as shown in Figure 10a. A constant pressure is applied to the block for 0.12 s as depicted in Figure 10.

**Figure 10.** Problem definition of beam with a flying block.

The frictionless contact condition (2) is imposed to allow for minimal sliding in the tangential direction after the occurrence of the impact. Therefore, a situation in which the configurations on the non-matching interface inevitably occur is considered. The material properties for both the structures are assumed to be those of the structural steel (*<sup>E</sup>* = <sup>2</sup> × <sup>10</sup>11, Poisson's ratio is 0.3, and density *<sup>ρ</sup>* = 7890). The mean squared displacement of a point "A" is used as an objective function (17) and the volume is limited to 50% for the topology optimization process (15).

The optimization was performed considering only the lower beam and not the upper block. A dynamic analysis is preferentially conducted to confirm the contact kinematics before performing topology optimization. Figure 11 illustrates the deformed shapes obtained by conducting the dynamic analysis, where it is observed that the impact occurs after 0.06 s, and then the frictionless contact occurs between two parts.

The topology optimization results based on the analysis are presented in Figure 12. In the numerical example, the optimization is performed for three cases: no contact (t = 0.04), impact (t = 0.08) and steady state (t = 0.12). The topological change in the design domain does not occur since there is no contact until 0.06 s, as shown in Figure 12a. The optimized layouts in the impact (t = 0.08) and steady state (t = 0.12) are derived as shown in Figure 12b,c, after 0.06 s to minimize the beam deflection. A large portion of the materials are distributed in the position where the beam is impacted when the design domain is in the impact state (t = 0.08), since the force is not completely transmitted to the beam, as shown in Figure 12b. However, when the analysis is performed over a period of 0.12 s, a shape in which the material is spread throughout the design domain is observed since the load is distributed throughout the beam, as shown in Figure 12c. The result shown in Figure 12c is a truss-type layout similar to the optimized results that can be derived under static analysis with the same boundary conditions.

**Figure 11.** Deformed configuration for beam with flying block.

D2SWLPXPOD\RXWIRUW

E2SWLPXPOD\RXWIRUW

F2SWLPXPOD\RXWIRUW

**Figure 12.** Optimized for beam with flying block using condensed localized mortar method.

As a result, it is possible to obtain the optimum topological layout to which the dynamic effect is applied in the impact state in which the non-matching interface inevitably occurred. In addition, after the impact state, it is confirmed that the identical layout with the quasi-static state can be obtained when the load by the frictionless contact condition is transmitted to the entire design domain.

It is confirmed through this example that the proposed approach can derive an acceptable result in the case of dynamic analysis where impact and frictionless contact are observed. Consequently, it is concluded that this approach can be used to perform topology optimization for various linear dynamic impact and frictionless contact problems.

#### **5. Conclusions**

This study presents an approach to apply the treatment of the non-matching interface to the dynamic response topology optimization for determining an optimum topological layout in the problem where the non-matching interface inevitably exists owing to the linear dynamic impact and frictionless contact conditions. First, the condensed localized mortar method is used to deal with linear dynamic impact and frictionless contact conditions which must occur in the non-matching interface during the analysis. The detailed explanation is in Section 2, and its usefulness is verified through the verification examples in Section 2.3. Second, the SIMP method, the most widely used topology optimization approach, is applied to treat the linear dynamic problems with the non-matching interface by using the condensed localized mortar method. The detailed explanation is in Section 3, and its effectiveness is verified through the verification examples in Section 3.2. Since the non-matching interface could be considered as the general structural analysis in the same framework by the features of the condensed localized mortar method, it can be easily applied to the well-known and well-established linear dynamic optimization framework. The effectiveness of the proposed approach is verified through a numerical example in Section 4. From Section 4, it is confirmed that physically reasonable results are obtained when the presented approach is applied in the dynamic response topology optimization in which the non-matching interface inevitably exist during the impact and frictionless contact conditions. Consequently, the problem of relative behavior including the impact and frictionless contact conditions along the adjacent parts could be adequately designed using the proposed method.

This study focuses on the linear dynamic impact and frictionless contact conditions to determine an optimum topological layout. In the future, the proposed optimization approach will be applied to more general and highly nonlinear non-matching interface problems, such as friction contact and multi-physics problems.

**Funding:** This work was supported by the Institute for Korea Spent Nuclear Fuel (iKSNF) and National Research Foundation of Korea (NRF) grant funded by the Korean government (Ministry of Science and ICT, MSIT) (No. 2021M2E1A1085229).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The author specially would like to thanks Sung-Kie Youn (Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology) for useful discussion and suggestion and constant academic supports.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviations/Nomenclature:**


#### **References**

