\jyear

2022

\nomail

[1]\fnmTsuyoshi \surNomura

\nomail\equalcont

Present address: \orgnameToyota Physical and Chemical Research Institute, \orgaddress41-1, Yokomichi, Nagakute, \cityAichi \postcode480-1192, \countryJapan

1]\orgnameToyota Central R&D Labs., Inc., \orgaddress41-1, Yokomichi, Nagakute, \cityAichi \postcode480-1192, \countryJapan

Topology optimization of locomoting soft bodies using material point method

\fnmYuki \surSato \fnmHiroki \surKobayashi \fnmChangyoung \surYuhn \fnmAtsushi \surKawamoto nomu2@mosk.tytlabs.co.jp \fnmNoboru \surKikuchi [

Abstract

Topology optimization methods have widely been used in various industries, owing to their potential for providing promising design candidates for mechanical devices. However, their applications are usually limited to the objects which do not move significantly due to the difficulty in computationally efficient handling of the contact and interactions among multiple structures or with boundaries by conventionally used simulation techniques. In the present study, we propose a topology optimization method for moving objects incorporating the material point method, which is often used to simulate the motion of objects in the field of computer graphics. Several numerical experiments demonstrate the effectiveness and the utility of the proposed method.

keywords:

Topology optimization, Material point method, Soft body, Soft robotics

1 Introduction

Topology optimization methods have successfully been used in many industries, such as automotive industries, owing to their potential for obtaining promising design candidates at the conceptual design stage. The key concept of topology optimization techniques is the replacement of structural optimization with a material distribution problem, which allows the topological changes in addition to changes in shapes during optimization procedures. After a pioneering work by Bendsøe and Kikuchi (1988), topology optimization techniques have been applied to various design problems, including multi-physics problems (Sigmund, 2001; Dilgen et al, 2018) and manufacturing-oriented design problems (Liu et al, 2018). Topology optimization techniques have also been applied to dynamical design problems, such as linkage mechanisms (Han et al, 2017; Han and Kim, 2021) and gear-linkage mechanisms (Yim et al, 2019). However, topology optimization methods have rarely been applied to design the structure of locomoting soft bodies, which involve large deformations of materials themselves. It would be due to the difficulty in handling the contact and interaction among multiple structures or with boundaries considering the geometrical non-linearity in a computationally efficient manner by conventionally used simulation techniques represented by the finite element method (FEM).

Handling locomoting soft bodies is essential in the field of soft robotics for developing deformable robots. This field has recently attracted much attention due to the relation to multidisciplinary fields, covering electronics, materials science, computer science, and biomechanics (Laschi and Cianchetti, 2014; Whitesides, 2018), which implies that the field has the potential for providing various applications. In the literature, several previous studies proposed controller optimization of locomoting soft robots (George Thuruthel et al, 2018; Kim et al, 2021; Hu et al, 2020; Bhatia et al, 2021). Additionally, Cheney et al (2014) performed design optimization of locomoting soft robots using the compositional pattern-producing network (CPPN). Bhatia et al (2021) also performed design optimization of locomoting soft robots based on genetic algorithms and Bayesian optimization in addition to CPPN. However, since the number of design variables this method can efficiently handle is relatively small, the degree of design freedom is limited. Therefore, there is room for improving the performance of soft robots by further optimizing their structures. Furthermore, this opens up a new possibility of finding a novel soft mechanism by computational design.

For topology optimization of locomoting soft bodies, numerical simulation techniques, which are computationally efficient in handling the dynamics of soft bodies, are required. One of the possible candidates is the material point methods (MPMs) (Sulsky et al, 1994; Jiang et al, 2016), which were originally introduced in computational mechanics and are now widely used to efficiently compute the motion of moving objects in the field of computer graphics (Stomakhin et al, 2013, 2014). MPMs are the hybrid Lagrangian–Eulerian techniques and have advantages over conventional Lagrangian or Eulerian methods. Compared with FEMs, which are conventionally used in topology optimization, the MPM can easily handle large deformation, collision, and contact, owing to the hybrid Lagrangian–Eulerian approach.

When using the FEMs, the large deformation simulations often cause numerical instabilities due to excessive distortion of meshes having low stiffness (Wang et al, 2014). On the other hand, MPMs easily simulate large deformation because they are mesh-free approaches, tracking the deformation using particles in the Lagrangian description. Furthermore, MPMs also simulate the collision and contact relatively easily by calculating the interaction on grids in the Eulerian description. Thus, the hybrid Lagrangian–Eulerian approach makes it easier to simulate large deformation, collision, and contact.

Recently, a topology optimization method incorporating MPMs has been proposed to optimize elastic materials with large static deformation (Li et al, 2021). Another technique to be required for topology optimization of locomoting soft bodies is the topology optimization for time-domain problems. Wang et al (2020) proposed space–time topology optimization where the time-dependent material distribution is optimized for additive manufacturing. Nomura et al (2007) and Yaji et al (2018) formulated the optimization problems in the time domain where the forward and sensitivity analyses were performed in a relevant time interval. In the present study, we propose a topology optimization method for locomoting soft bodies incorporating MPMs, formulating the optimization problem in the time domain. The proposed method incorporates a material representation scheme for topology optimization with MPMs to optimize locomoting soft bodies as material distribution, with a formulation of a density filtering in the framework of MPMs for obtaining smoothed structures. The contributions of the present study are the followings:

•

We incorporated the MPM, originally introduced in computational mechanics but now widely used in computer graphics, with structural optimization methods, especially topology optimization methods. This incorporation enables topology optimization of locomoting soft bodies, which might otherwise be challenging for conventional schemes.
•

We introduced topology optimization, i.e., material distribution optimization, into the design optimization of locomoting soft bodies. This makes the degree of the design freedom much higher compared to the existing approaches for designing soft bodies.
•

We fully exploited the framework of DiffTaichi (Hu et al, 2020), especially for the sensitivity calculation of many design variables. The utilization of the automatic differentiation (AD) is essential for dynamic design problems, which involve complicated situations such as contacts.

2 Formulations

2.1 Governing equation

Firstly, the conservation law of the mass and momentum reads

	$\displaystyle\dfrac{D\rho}{Dt}+\rho\nabla\cdot\boldsymbol{v}=0,$		(1)
	$\displaystyle\rho\dfrac{D\boldsymbol{v}}{Dt}=\nabla\cdot\boldsymbol{\sigma}+\rho\boldsymbol{g},$		(2)

where $\rho$ is the density, $t$ is the time, $\boldsymbol{v}$ is the velocity, $\boldsymbol{\sigma}$ is the Cauchy stress tensor, and $\boldsymbol{g}$ is the gravity acceleration. The operator $D/Dt$ represents the material derivative. Let $\boldsymbol{F}$ denote the deformation gradient, then the Cauchy stress is given as

\boldsymbol{\sigma}=\frac{1}{J}\dfrac{\partial\varPsi}{\partial\boldsymbol{F}}\boldsymbol{F}^{\top},

(3)

where $J=\mathrm{det}(\boldsymbol{F})$ and $\varPsi$ is the potential energy. Assuming the neo-Hookean material constitution, the potential energy $\varPsi$ is given as follows:

\displaystyle\varPsi(\boldsymbol{F}):=\dfrac{\mu}{2}\left(\mathrm{tr}(\boldsymbol{F}^{\top}\boldsymbol{F})-d\right)-\mu\log(J)+\dfrac{\lambda}{2}\log^{2}(J),

(4)

where $d$ is the spatial dimensions. These governing equations are solved using the material point method (MPM). Specifically, in the present study, we use the moving least square material point method (MLS-MPM) (Hu et al, 2018), which incorporates the moving least square discretization of the governing equation with the MPM.

2.2 Material representation in material point method

We now describe the material representation in the material point method. Since the governing equations include a self-weight load, i.e. the gravity acceleration, certain treatments in material interpolation schemes are required for dealing with it. Several previous studies (Bruyneel and Duysinx, 2005; Kumar, 2022) proposed novel interpolation schemes of the mass density so that numerical instabilities at low densities could be avoided and the non-monotonous behaviors of objective functions could be reduced. In the present study, since the objective function is originally non-monotonous, different from structural compliance, the SIMP interpolation does not contribute to penalizing the intermediate values of design variables, i.e., grayscales. Thus, we used a linear interpolation scheme both for stiffness and mass density, which avoided the ratio of the gravity force to stiffness from diverging and made it easier to handle self-weight problems. Let $\gamma_{p}$ for $p=1,\ldots,N_{p}$ denote the fictitious material density of the $p$ -th particle in the design domain where $N_{p}$ is the number of particles in the domain. Then, the density and the Lamé’s constants of the $p$ -th particle are linearly interpolated using $\gamma_{p}$ as follows:

$\displaystyle\rho_{p}$	$\displaystyle=\rho\left((1-\varepsilon)\gamma_{p}+\varepsilon\right),$	(5)
$\displaystyle\lambda_{p}$	$\displaystyle=\lambda\left((1-\varepsilon)\gamma_{p}+\varepsilon\right),$	(6)
$\displaystyle\mu_{p}$	$\displaystyle=\mu\left((1-\varepsilon)\gamma_{p}+\varepsilon\right),$	(7)

where $\lambda$ and $\mu$ are the Lamé’s constants of the material, and $\varepsilon$ is a small constant, which is introduced to avoid the singularity.

Refer to caption — Figure 1: Schematic of the density filtering. In the particle-to-grid process, design variables of particles in the red region are transferred to the $(i,j)$ grid (blue node). In the grid operation, smoothed grid design variables are calculated as the convolution of the grid design variables. In the grid-to-particle process, smoothed grid design variables at $(i+k,j+l)$ grids for $k,l\in\{-1,0,1\}$ ( $3\times 3$ blue nodes) are transferred to smoothed design variables of particles in the blue region.

2.3 Heaviside projection with density filtering

To ensure the smoothness of the structures, topology optimization methods usually employ a certain filtering technique. In the present study, we introduce a filtering technique in the scheme of the material point method. Let $\phi_{p}$ denote the auxiliary design variables assigned to the $p$ -th particle to be designed. These design variables are spatially filtered as shown in Fig. 1 and projected onto the fictitious particle densities $\gamma_{p}$ , as follows:

1.

Particle to grid: Transfer the particle design variables $\phi_{p}$ to grid design variables using the same framework as the MLS-MPM. The transferred grid design variables are then normalized by the sum of the weights of each particle to each grid.

Grid operation: Take the weighted sum of the grid design variables in the neighboring grid points. For two-dimensional problems, for example, let $\phi_{i,j}$ and $\tilde{\phi}_{i,j}$ respectively denote the grid design variables and smoothed grid design variables on the $(i,j)$ grid. Then the smoothed grid design variables are calculated as the convolution of the grid design variables, as follows:

$\displaystyle\tilde{\phi}_{i,j}=$	$\displaystyle w_{0}\phi_{i,j}+w_{1}\phi_{i-1,j}+w_{1}\phi_{i+1,j}$
	$\displaystyle+w_{1}\phi_{i,j-1}+w_{1}\phi_{i,j+1}+w_{2}\phi_{i-1,j-1}$
	$\displaystyle+w_{2}\phi_{i+1,j-1}+w_{2}\phi_{i-1,j+1}+w_{2}\phi_{i+1,j+1},$	(8)

where $w_{0},w_{1},w_{2}$ are weighting coefficients satisfying $w_{0}+4w_{1}+4w_{2}=1$ . These parameters $w_{0}$ , $w_{1}$ , and $w_{2}$ control the degree of smoothness of filtered variables.

3.

Grid to particle: Transfer the smoothed grid design variables to the smoothed particle design variables $\tilde{\phi}_{p}$ using the same framework as the MLS-MPM.
4.

Heaviside projection: Project the smoothed particle design variables $\tilde{\phi}_{p}$ onto the fictitious particle densities using the smoothed Heaviside function as follows:

$\gamma_{p}=\dfrac{1}{2}\left(\dfrac{\tanh{(\beta\tilde{\phi}_{p})}}{\tanh{(\beta})}+1\right),$ (9)

where $\beta$ is a parameter to control the curvature of the smoothed Heaviside function.

2.4 Problem statement

In the present study, we focus on the problem of designing soft robots that move toward the designated direction as far as possible. Consider that a soft robot to be designed consists of the fixed design domain, the actuation domain, and the non-design domain. In the fixed design domain, the material distribution is optimized, while in the other domains, the structure is unchanged during the optimization process. In the actuation domain, the prescribed stress $\boldsymbol{\sigma}_{\mathrm{act}}$ is added to the Cauchy stress in Eq. (3). The prescribed stress $\boldsymbol{\sigma}_{\mathrm{act}}$ is given as

\boldsymbol{\sigma}_{\mathrm{act}}(t)=\sum_{k=1}^{N_{a}}\left(a_{k}\sin(\omega_{k}t+\pi\theta_{k})+b_{k}\right)\boldsymbol{F}(t)\boldsymbol{S}_{k}\boldsymbol{F}(t)^{\top},

(10)

where $\boldsymbol{S}_{k}$ is the second Piola-–Kirchhoff stress caused by the $k$ -th actuator, $N_{a}$ is the number of actuators, $a_{k}$ is the amplitude, $\theta_{k}$ is the phase, $b_{k}$ is the bias, and $\omega_{k}$ is angular frequency.

The objective is to maximize the distance that an object travels, which is measured by the mass center in the designated direction $\boldsymbol{e}$ at the end of the prescribed time interval as follows:

\displaystyle\mathcal{L}:=\dfrac{\int_{\Omega_{N}(T)}\boldsymbol{x}\cdot\boldsymbol{e}~{}d\boldsymbol{x}}{\int_{\Omega_{N}(T)}d\boldsymbol{x}},

(11)

where $\Omega_{N}(T)$ represents the union of the non-design and actuation domains at predetermined time $T$ . Here, we define the objective function within these domains so that the domain of integration should be unchanged during the optimization procedure. The design variables are $\boldsymbol{\phi}=\left[\phi_{1},\ldots,\phi_{N_{p}}\right]^{\top}$ for material configurations and the input parameters to actuators $\boldsymbol{a}=\left[a_{1},\ldots,a_{N_{a}}\right]^{\top}$ , $\boldsymbol{\theta}=\left[\theta_{1},\ldots,\theta_{N_{a}}\right]^{\top}$ , and $\boldsymbol{b}=\left[b_{1},\ldots,b_{N_{a}}\right]^{\top}$ . All design variables are bounded within $[-1,1]$ . The optimization problem is then formulated in the minimization form as


$\displaystyle\min_{\boldsymbol{\phi},\boldsymbol{a},\boldsymbol{\theta},\boldsymbol{b}}$	$\displaystyle-\mathcal{L},$	(12a)
subject to	$\displaystyle\boldsymbol{\phi}\in[-1,1]^{N_{p}},$	(12b)
$\displaystyle\boldsymbol{a},\boldsymbol{\theta},\boldsymbol{b}\in[-1,1]^{N_{a}},$		(12c)

in which Eqs. (1) and (2) are nested.

3 Numerical experiments

In the following, we provide two numerical experiments. The proposed method was implemented by Python and specifically, the forward and backward analysis using the MPM was implemented by Taichi (Hu et al, 2019, 2020), which can be called in Python environment. Taichi is an open-source programming language for high-performance numerical computation, and provides a scheme for automatic differentiation which we used for sensitivity calculations. Lamé’s constants were respectively set to $\lambda=10$ and $\mu=10$ , and the small constant $\varepsilon$ was set to $10^{-3}$ . The density was set to $\rho=1$ . The number of actuators $N_{a}$ was set to $1$ , and the angular frequency of the actuators $\omega_{1}$ was set to $30$ . The prescribed second Piola-Kirchhoff stress was given as $S_{1}=2I$ where $I$ is the identity matrix. The weighting coefficients for density filtering were set to $w_{0}=0.64$ , $w_{1}=0.08$ , and $w_{2}=0.01$ . These parameters $w_{0}$ , $w_{1}$ , and $w_{2}$ control the degree of smoothness of filtered variables and were determined empirically through numerical experiments. The number of grids in the MPM was $128\times 128$ , and the particles were initially aligned in a grid-like manner with a doubled resolution to the grids. The forward simulation was performed $5000$ time steps with the time interval $2.0\times 10^{-4}$ s. The optimization was performed using the method of moving asymptotes (Svanberg, 1987). The parameter $\beta$ in the smoothed Heaviside function was initially set to $0.5$ and was doubled when the optimization process converged, until the parameter reached $16$ . The optimization process continued while the change in the moving average of the objective function values over $10$ iterations was larger than $2.0\times 10^{-4}$ .

3.1 Walker design

The first problem setting for designing a walking robot is illustrated in Fig. 2. The gravity acceleration $\boldsymbol{g}$ is set to $-10$ along the $y$ -axis. The optimization result is shown in Fig. 3, which exhibits that the obtained structure can walk toward the designated direction, i.e., along the $x$ -axis. The input parameters of the actuator were optimized to $a=1.00$ , $\theta=-0.330$ , and $b=-0.982$ . In the optimized structure, legs-like substructures were formed, which took a step toward the designated right direction as shown in Fig. 3 (b) and (c). To verify the optimized result, we further performed forward simulations of the walker’s body exactly represented by particles placed where $\rho\geq 0.3$ , which was empirically set. The post-processed results still exhibit better performance compared to the reference designs shown in Figs. 4. The objective function values were listed in Table 1.

Table 1: Objective function values of the walker.

Condition	Objective function value
Optimized	-0.362
Post-processed	-0.357
Reference design	-0.208

3.2 Crawler design

The second problem setting for designing a crawling robot is illustrated in Fig. 5. In this problem, no gravity is applied. Moreover, the grid design variables in the lower half are mapped to those in the upper half just after Step 2 in Sec. 2.3 to keep the design symmetry across the horizontal mid-axis. Figure 6 shows the optimized result of the crawler design. The input parameters of the actuator were optimized to $a=0.964$ , $\theta=-0.174$ , and $b=-0.447$ . As shown in this figure, the optimized structure moves forward by alternatively kicking the walls with the front and rear legs. As a verification of the results, we performed forward simulations in which the crawler’s body was exactly represented by particles placed where $\rho\geq 0.3$ . The post-processed results exhibit better performance compared to the reference designs shown in Fig. 7 as listed in Table 2.

Table 2: Objective function values of the crawler.

Condition	Objective function value
Optimized	-0.451
Post-processed	-0.401
Reference design	-0.379

4 Conclusion

This paper presented a topology optimization for locomoting soft bodies incorporating the material point method. We proposed a material representation scheme in MPMs and constructed a density filtering technique in the framework of MPMs, which enables us to use MPMs in topology optimization. The effectiveness of the proposed method was confirmed through numerical experiments. On the other hand, the present study has limitations on manufacturability due to the material property setting. Since the dynamics we focus on are geometrically non-linear, the material property and geometry must be set based on the real-world environment. Albeit such requirements, the present study used numerically convenient settings for simplicity. We hope to address manufacturability issues in our future research. We believe that the proposed method will broaden the field of applications of topology optimization for the dynamics of soft bodies, which is related to large deformation, collision, and contact.

Statements and Declarations

\bmhead

Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest. \bmheadReplication of results The source code is unavailable due to institutional constraints. However, further algorithm details are available upon request to the authors.

Acknowledgement

The manuscript has no funding for this paper or is not applicable.

References

\bibcommenthead
Bendsøe and Kikuchi (1988) Bendsøe MP, Kikuchi N (1988) Generating optimal topologies in structural design using a homogenization method. Computer methods in applied mechanics and engineering 71(2):197–224
Bhatia et al (2021) Bhatia J, Jackson H, Tian Y, Xu J, Matusik W (2021) Evolution gym: A large-scale benchmark for evolving soft robots. Advances in Neural Information Processing Systems 34:2201–2214
Bruyneel and Duysinx (2005) Bruyneel M, Duysinx P (2005) Note on topology optimization of continuum structures including self-weight. Structural and Multidisciplinary Optimization 29(4):245–256
Cheney et al (2014) Cheney N, MacCurdy R, Clune J, Lipson H (2014) Unshackling evolution: evolving soft robots with multiple materials and a powerful generative encoding. ACM SIGEVOlution 7(1):11–23
Dilgen et al (2018) Dilgen SB, Dilgen CB, Fuhrman DR, Sigmund O, Lazarov BS (2018) Density based topology optimization of turbulent flow heat transfer systems. Structural and Multidisciplinary Optimization 57(5):1905–1918
George Thuruthel et al (2018) George Thuruthel T, Ansari Y, Falotico E, Laschi C (2018) Control strategies for soft robotic manipulators: a survey. Soft robotics 5(2):149–163
Han and Kim (2021) Han SM, Kim YY (2021) Topology optimization of linkage mechanisms simultaneously considering both kinematic and compliance characteristics. Journal of Mechanical Design 143(6)
Han et al (2017) Han SM, In Kim S, Kim YY (2017) Topology optimization of planar linkage mechanisms for path generation without prescribed timing. Structural and Multidisciplinary Optimization 56(3):501–517
Hu et al (2018) Hu Y, Fang Y, Ge Z, Qu Z, Zhu Y, Pradhana A, Jiang C (2018) A moving least squares material point method with displacement discontinuity and two-way rigid body coupling. ACM Transactions on Graphics 37(4):150
Hu et al (2019) Hu Y, Li TM, Anderson L, Ragan-Kelley J, Durand F (2019) Taichi: a language for high-performance computation on spatially sparse data structures. ACM Transactions on Graphics (TOG) 38(6):201
Hu et al (2020) Hu Y, Anderson L, Li TM, Sun Q, Carr N, Ragan-Kelley J, Durand F (2020) DiffTaichi: differentiable programming for physical simulation. ICLR
Jiang et al (2016) Jiang C, Schroeder C, Teran J, Stomakhin A, Selle A (2016) The material point method for simulating continuum materials. In: ACM SIGGRAPH 2016 Courses. p 1–52
Kim et al (2021) Kim D, Kim SH, Kim T, Kang BB, Lee M, Park W, Ku S, Kim D, Kwon J, Lee H, et al (2021) Review of machine learning methods in soft robotics. PLoS One 16(2):e0246,102
Kumar (2022) Kumar P (2022) Topology optimization of stiff structures under self-weight for given volume using a smooth heaviside function. Structural and Multidisciplinary Optimization 65(4):1–17
Laschi and Cianchetti (2014) Laschi C, Cianchetti M (2014) Soft robotics: new perspectives for robot bodyware and control. Frontiers in bioengineering and biotechnology p 3
Li et al (2021) Li Y, Li X, Li M, Zhu Y, Zhu B, Jiang C (2021) Lagrangian–Eulerian multidensity topology optimization with the material point method. International Journal for Numerical Methods in Engineering 122(14):3400–3424
Liu et al (2018) Liu J, Gaynor AT, Chen S, Kang Z, Suresh K, Takezawa A, Li L, Kato J, Tang J, Wang CC, et al (2018) Current and future trends in topology optimization for additive manufacturing. Structural and multidisciplinary optimization 57(6):2457–2483
Nomura et al (2007) Nomura T, Sato K, Taguchi K, Kashiwa T, Nishiwaki S (2007) Structural topology optimization for the design of broadband dielectric resonator antennas using the finite difference time domain technique. International Journal for Numerical Methods in Engineering 71(11):1261–1296
Sigmund (2001) Sigmund O (2001) Design of multiphysics actuators using topology optimization–part I: one-material structures. Computer methods in applied mechanics and engineering 190(49-50):6577–6604
Stomakhin et al (2013) Stomakhin A, Schroeder C, Chai L, Teran J, Selle A (2013) A material point method for snow simulation. ACM Transactions on Graphics (TOG) 32(4):1–10
Stomakhin et al (2014) Stomakhin A, Schroeder C, Jiang C, Chai L, Teran J, Selle A (2014) Augmented MPM for phase-change and varied materials. ACM Transactions on Graphics (TOG) 33(4):1–11
Sulsky et al (1994) Sulsky D, Chen Z, Schreyer HL (1994) A particle method for history-dependent materials. Computer methods in applied mechanics and engineering 118(1-2):179–196
Svanberg (1987) Svanberg K (1987) The method of moving asymptotes–a new method for structural optimization. International journal for numerical methods in engineering 24(2):359–373
Wang et al (2014) Wang F, Lazarov BS, Sigmund O, Jensen JS (2014) Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems. Computer Methods in Applied Mechanics and Engineering 276:453–472
Wang et al (2020) Wang W, Munro D, Wang CC, van Keulen F, Wu J (2020) Space-time topology optimization for additive manufacturing. Structural and Multidisciplinary Optimization 61(1):1–18
Whitesides (2018) Whitesides GM (2018) Soft robotics. Angewandte Chemie International Edition 57(16):4258–4273
Yaji et al (2018) Yaji K, Ogino M, Chen C, Fujita K (2018) Large-scale topology optimization incorporating local-in-time adjoint-based method for unsteady thermal-fluid problem. Structural and Multidisciplinary Optimization 58(2):817–822
Yim et al (2019) Yim NH, Kang SW, Kim YY (2019) Topology optimization of planar gear-linkage mechanisms. Journal of Mechanical Design 141(3):032,301