A Two-Player Game for Multi-Scale Topology Optimization of Static and Dynamic Compliances of Triply Periodic Minimal Surface-Based Lattice Structures

Abstract

In this study, a novel non-cooperative two-player game for minimizing static (Player 1) and dynamic (Player 2) compliances is introduced, implemented, and demonstrated using a multi-scale topology optimization framework for triply periodic minimal surface (TPMS)-based lattice structures. Player 1 determines the optimal macro-layout by minimizing the static compliance based on a micro-layout provided by Player 2. Conversely, player 2 identifies the optimal micro-layout (grading of the TPMS-based lattice structure) by minimizing the dynamic compliance given a macro-layout from Player 1. The multi-scale topology optimization formulations are derived using two density variables in each finite element. The first variable is the standard density, which dictates whether the finite element is void or contains the graded lattice structure and is governed by the rational approximation of material properties (RAMP) model. The second density variable represents the local relative density of the TPMS-based lattice structure, determining the effective orthotropic elastic properties of the finite element. The multi-scale game is implemented for three-dimensional problems, and solved using a Gauss–Seidel algorithm with sequential linear programming. It is numerically demonstrated for several benchmarks that the proposed multi-scale game generates equilibrium designs with strong performance for both static and harmonic load cases, effectively avoiding resonance at harmonic load frequencies. Validation is achieved through modal analyses of finite element models of the optimal designs.

1. Introduction

Additive manufacturing enables the production of complex, cellular-based, lightweight multifunctional designs that exhibit both strength and stiffness [1,2]. These designs are similar to lightweight structures appearing in nature, such as bird skeletons [3] and butterfly wings [4]. A particularly suitable class of self-supporting cellular structures for this purpose is the triply periodic minimal surface (TPMS)-based lattice structures, developed and described by Schoen [5], with the Gyroid structure being the most well known.

TPMS-based structures are configured as implicit surface-based geometries using Boolean operators, resulting in either frame-based or shell-based structures, which are straightforward to grade [6]. Recently, Strömberg [7] proposed a multi-scale topology optimization framework for determining optimal macro-layouts with optimal grading of TPMS-based lattice structures possessing transversely isotropic elastic bulk properties. In an earlier study, Strömberg [8] implemented frame- and shell-based Gyroid, G-prime, and Schwarz-D structures by establishing elastic properties as a function of relative lattice density using numerical homogenization of representative volume elements of the TPMS-based lattice structures. In this study, those previous works [7,8] are further developed by including harmonic excitations. A recent review on multi-scale topology optimization can be found in [9].

When the frequency of a harmonic load approaches any of the natural frequencies (eigenfrequencies) of a component, resonance occurs. This phenomenon is illustrated in Figure 1 for a compliance optimal pin-jointed beam, where a critical natural frequency at 1604 Hz is identified using modal analysis. Consequently, if an external harmonic force has an excitation frequency near 1604 Hz, resonance will occur, rendering the optimal compliance design suboptimal for such harmonic loads. Therefore, in the design of lightweight components, natural frequencies close to the driving frequency of the external load must be avoided. One design approach to achieve this involves utilizing TPMS-based lattice structures and grading them to tune the natural frequencies, thereby avoiding resonance [10]. This concept is explored in the following paper by proposing a novel multi-scale two-player topology optimization game, by adopting the multi-scale topology optimization framework for TPMS-based lattice structures recently proposed by Strömberg [7].

The concept of minimizing dynamic compliance in topology optimization can be traced back to the work by Ma et al. [11]. Other notable early contributions include papers by Jog [12], and Olhoff and Du [13]. Dynamic compliance is defined as $|{\mathit{F}}^T\mathit{u}|$, where $\mathit{F}$ is the amplitude of the external harmonic force, and $\mathit{u}$ represents the magnitudes of the displacement response obtained by solving

$$\left(-{\Omega }^2\mathit{M}+\mathit{K}\right)\mathit{u}=\mathit{F},$$

where $\Omega$ is the harmonic excitation angular frequency, and $\mathit{M}$ and $\mathit{K}$ are the mass and stiffness matrices, respectively. Thus, by minimizing the dynamic compliance, resonances at critical eigenmodes can be avoided by tuning $\mathit{M}$ and $\mathit{K}$. In the proposed two-player game, this kind of tuning will be accomplished by finding optimal gradings (micro-layouts) of the TPMS-based lattice structures by minimizing the square of the dynamic compliance.

A potential drawback of the dynamic compliance approach is that the design may disconnect from the supports, resulting in rigid body modes becoming the preferred design solutions. This leads to the static compliance approaching infinity [14]. This draw-back is also discussed in the textbook by Bendsœ and Sigmund [15]. This issue can be mitigated by including a constraint on static compliance in the topology optimization formulation. In the proposed two-player game, the development of rigid body modes is instead prevented by Player 1, who finds optimal macro-layouts of TPMS-based lattices by minimizing the static compliance.

It can be shown that the dynamic compliance formulation forces the closest eigenfrequency away from the excitation frequency [15]. Thus, the dynamic compliance formulation should not be confused with the maximum eigenvalue formulations, as discussed in the introduction by Tsai and Cheng [16] and in [15]. Further discussions on the selection of objective functions for minimizing the structural dynamic response were presented by Niu et al. [17]. Silva et al. [18] recently presented a critical analysis on this topic for single-material structures. The multi-scale approach adopted in this work can also be viewed as a multi-material approach. Liu et al. [19] presented recent work on bi-material structures for maximizing band-gap. Long et al. [20] performed topology optimization of transient dynamic problems using the second-order Arnoldi reduction scheme. Zargham et al. [21] conducted a review on topology optimization of structural vibration problems, while Simsek et al. [22] studied the structural dynamic performance of Gyroid structures.

A non-cooperative two-player game consists of two players, denoted as player $P_1$ and player $P_2$, each controlling their own strategy with variables collected in ${\mathit{x}}_i$ ($i=1,2$). Each players aims to find a strategy that maximize its own payoff $f_i=f_i({\mathit{x}}_1,{\mathit{x}}_2)$ without collaborating with the other player. An optimal strategy is denoted ${\mathit{x}}_i^{*}$. A pair of optimal strategies $({\mathit{x}}_1^{*},{\mathit{x}}_2^{*})$ is termed a Nash equilibrium if the following conditions hold:

$$\begin{array}{c}{P}_1:f_1({\mathit{x}}_1,{\mathit{x}}_2^{*})\le f_1({\mathit{x}}_1^{*},{\mathit{x}}_2^{*})\hspace{1em}\forall {\mathit{x}}_1,\hfill \\ P_2:f_2({\mathit{x}}_1^{*},{\mathit{x}}_2)\le f_2({\mathit{x}}_1^{*},{\mathit{x}}_2^{*})\hspace{1em}\forall {\mathit{x}}_2.\hfill \end{array}$$

The two-player game approach for topology optimization utilized in this work draws significant inspiration from the early work by Habbal et al. [23]. Subsequently, Holmberg et al. [24] employed this approach for robust topology optimization under uncertainty in loading conditions. Thore et al. [25] further investigated structural optimization under uncertainty using similar two-player game approaches. Most recently, Thore et al. [26] established a similar mathematical game for topology optimization of cooling systems. A review of multi-objective optimization and its applications can be found in [27].

The recent study presented in [28] reviews the use of additive manufacturing, topology optimization, and lattice structures in lightweight design, highlighting various industrial applications. In the aerospace and automotive sectors, where lightweight components are critical, fatigue often serves as the primary design constraint. The work in [29] specifically examines design optimization and additive manufacturing of a fatigue-critical aerospace component, while fatigue failure predictions in lattice structures are explored in [30]. When a harmonic excitation approaches an eigenfrequency, the risk of fatigue failure becomes significantly higher, making it crucial to avoid resonance by tuning the eigenfrequencies away from the excitation frequency. The two-player game proposed in this paper addresses this issue by grading the TPMS-based lattice structures, effectively mitigating the risk of resonance. The numerical implementation of the proposed game using a Gauss–Seidel algorithm with sequential linear programming (SLP) seems most promising for real-world engineering applications, as the computational cost is only marginally higher—approximately double—than solving a standard compliance topology optimization problem. The computational efficiency of the proposed method is demonstrated through the solution of a large-scale benchmark problem [31] presented at the end of the paper.

The outline of the paper is as follows: in the next section, the governing equations required for the multi-scale topology optimization formulation are presented. Section 3 formulates the strategies for the multi-scale players as two separate topology optimization problems. Section 4 presents the treatment of the multi-scale game using a Gauss–Seidel algorithm with SLP. Section 5 solves numerical benchmarks in 3D using the proposed non-cooperative two-player game. Finally, some concluding remarks are provided.

2. Governing Equations

In this section, the governing equations presented in [7] are reviewed and modified in order to include inertia. The modified governing equations are then adopted in the next section for defining the topology optimization problems used to set up the new two-player game.

A non-homogenous linear orthotropic elastic body with fixed displacements and prescribed external forces, which is discretized with linear finite elements, is considered. Two density variables, ${\gamma }_e$ and ${\rho }_e$, are introduced for each finite element e; ${\gamma }_e$ is the local relative density of lattice in the finite element, and ${\rho }_e$ is a standard macro-density variable defining if the finite element should be treated as a void or filled with lattice structure.

The effective elastic properties of each finite element e are obtained by numerical homogenization of representative volume elements of the TPMS-based lattice structure of interest such that the effective elasticity matrix in Voigt notation is given by

$${\mathit{C}}_e={\mathit{C}}_e\left({\gamma }_e\right)=\left[\begin{array}{cccccc}{f}_{11}{c}_{11}& f_{21}{c}_{21}& f_{31}{c}_{31}& 0& 0& 0\\ f_{21}{c}_{21}& f_{22}{c}_{11}& f_{32}{c}_{31}& 0& 0& 0\\ f_{31}{c}_{31}& f_{32}{c}_{31}& f_{33}{c}_{33}& 0& 0& 0\\ 0& 0& 0& f_{44}{c}_{44}& 0& 0\\ 0& 0& 0& 0& f_{55}{c}_{44}& 0\\ 0& 0& 0& 0& 0& f_{66}{c}_{66}\end{array}\right],$$

(1)

where $f_{ij}=f_{ij}\left({\gamma }_e\right)$ are material interpolation (regression) laws as a function of the relative lattice density $0\le l_b\le {\gamma }_e\le u_b\le 1$ of element e, and where $l_b$ and $u_b$ are prescribed lower and upper limits on ${\gamma }_e$, respectively [8]. The values of these limits are set to promote 3D printing of TPMS-based lattices with good quality. In the numerical examples, we let $l_b=0.2$ and $u_b=0.6$, which is an example of one such conservative setting.

Furthermore, $c_{ij}$ in (1) are the elastic properties of the bulk material of the lattice structure, which are assumed to be governed by transversely isotropic elasticity, i.e., ${\mathit{C}}_e{\left(1\right)}^{-1}=$

$$\left[\begin{array}{cccccc}1/E& -\nu /E& -{\nu }_{13}/E& 0& 0& 0\\ -\nu /E& 1/E& -{\nu }_{13}/E& 0& 0& 0\\ -{\nu }_{13}/E& -{\nu }_{13}/E& 1/E_{33}& 0& 0& 0\\ 0& 0& 0& 1/G_{23}& 0& 0\\ 0& 0& 0& 0& 1/G_{23}& 0\\ 0& 0& 0& 0& 0& 2(1+\nu )/E\end{array}\right],$$

(2)

where E and $E_{33}$ are Young’s moduli, $\nu$ and ${\nu }_{13}$ are Poisson’s ratios, and $G_{23}$ is the out-of-plane shear moduli.

The corresponding local effective finite element stiffness matrix for element e is denoted by

$${\mathit{k}}_e={\mathit{k}}_e\left({\gamma }_e\right).$$

(3)

In the global assembly procedure, the RAMP model, developed by Stolpe and Svanberg [32], is utilized in order to define if the finite element e should be considered to be a void or filled with lattice structure. Thus, the global stiffness matrix is generated by

$$\mathit{K}=\mathit{K}(\mathit{\rho },\mathit{\gamma })=\bigcap _e{\displaystyle \frac{{\rho }_e}{1+n(1-{\rho }_e)}}{\mathit{k}}_e\left({\gamma }_e\right),$$

(4)

where $\mathit{\rho }=\left\{{\rho }_e\right\}$, $\gamma =\left\{{\gamma }_e\right\}$, ⋂ represents an assembly operator, n is the RAMP factor, and $ϵ\le {\rho }_e\le 1$ is the relative macro-density of lattice structure for element e. $ϵ$ is a small number in order to avoid singular stiffness matrices.

The global mass matrix is given by

$$\mathit{M}=\mathit{M}(\mathit{\rho },\mathit{\gamma })=\bigcap _e{\displaystyle \frac{{\rho }_e}{1+exp(-{\lambda }_m({\rho }_e-{\widehat{\rho }}_m))}}{\gamma }_e{\mathit{m}}_e,$$

(5)

where ${\mathit{m}}_e$ is local element mass matrix for an element e with ${\rho }_e={\gamma }_e=1$,

$$f_m=f_m\left({\rho }_e\right)={\displaystyle \frac{{\rho }_e}{1+exp(-{\lambda }_m({\rho }_e-{\widehat{\rho }}_m))}}$$

(6)

appearing in (5) is a Sigmoid-like function enforcing the local mass towards zeros for ${\rho }_e<{\widehat{\rho }}_m$. In this work, ${\lambda }_m=50$ and ${\widehat{\rho }}_m=0.2$. This setting for $f_m$ is plotted in Figure 2. Thus, for low densities, the mass converges quicker to zero than the total volume of bulk material presented below in (7). This approach prevents spurious modes to be developed in the void region.

The total volume of bulk material generated in the assembly procedure is given by

$$V^{\mathrm{bulk}}=V^{\mathrm{bulk}}(\mathit{\rho },\mathit{\gamma })=\sum _e{\rho }_e{\gamma }_e{V}_e,$$

(7)

where $V_e$ represents the total volume of each element e when ${\rho }_e={\gamma }_e=1$. By using (7), the following different volume measures can be identified:

$$\begin{array}{c}{V}^{\mathrm{design}}=V^{\mathrm{bulk}}(\mathbf{1},\mathbf{1}),\\ V^{\mathrm{macro}}=V^{\mathrm{macro}}\left(\mathit{\rho }\right)=V^{\mathrm{bulk}}(\mathit{\rho },\mathbf{1}),\\ V^{\mathrm{lattice}}=V^{\mathrm{lattice}}\left(\mathit{\gamma }\right)=V^{\mathrm{bulk}}(\mathbf{1},\mathit{\gamma }),\end{array}$$

(8)

where $V^{\mathrm{design}}$ is the total volume of the design domain, $V^{\mathrm{macro}}$ is the total volume of the macro-layout filled with a lattice structure, $V^{\mathrm{lattice}}$ is the total volume of the lattice in the design domain, and $\mathbf{1}$ is a vector of ones. The volume of lattice $V_{\mathrm{macro}}^{\mathrm{lattice}}$ in the macro-layout of lattice structure $V^{\mathrm{macro}}$ is given by

$$V_{\mathrm{macro}}^{\mathrm{lattice}}=V^{\mathrm{lattice}}-(V^{\mathrm{design}}-V^{\mathrm{macro}})l_b,$$

(9)

which will converge towards the volume of bulk material $V^{\mathrm{bulk}}$, because ${\rho }_e=ϵ$ in the void region.

3. The Multi-Scale Players

In this section, the two multi-scale players with their strategies for the non-cooperative game presented in the next section are defined. The first player minimizes the static compliance $c_s$ by finding the optimal macro-layout ${\mathit{\rho }}^{*}$, and the second player minimizes half of the square of the dynamic compliance, denoted by $c_d$, by finding the optimal local grading ${\mathit{\gamma }}^{*}$ of the TPMS-based lattice structure. Taking the square of the dynamic compliance, instead of simply using the dynamic compliance $|{\mathit{F}}^T\mathit{u}|$ in the objective function, is a standard approach for obtaining a smooth objective function [14]. Furthermore, the square or the absolute value of the dynamic compliance is used because the dynamic stiffness might become negative definite for frequencies higher than the fundamental frequency [13].

Thus, the two players’ strategies, $P_1$ and $P_2$, are defined by the following multi-scale topology optimization problems:

$$P_1\left\{\begin{array}{c}\mathrm{Given} {\mathit{\gamma }}^{*}:\hfill \\ {\displaystyle \underset{(\mathit{\rho },\mathit{d})}{min}{c}_s={\mathit{F}}^T\mathit{d}}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\mathit{K}(\mathit{\rho },{\mathit{\gamma }}^{*})\mathit{d}=\mathit{F},\hfill \\ V^{\mathrm{macro}}\left(\mathit{\rho }\right)\le {\widehat{V}}^{\mathrm{macro}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \end{array}\right.\hfill \end{array}\right.$$

(10)

$$P_2\left\{\begin{array}{c}\mathrm{Given} {\mathit{\rho }}^{*}:\hfill \\ {\displaystyle \underset{(\mathit{\gamma },\mathit{u})}{min}{c}_d={\displaystyle \frac{1}{2}}{\left({\mathit{F}}^T\mathit{u}\right)}^2}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\left(-{\Omega }^2\mathit{M}({\mathit{\rho }}^{*},\mathit{\gamma })+\mathit{K}({\mathit{\rho }}^{*},\mathit{\gamma })\right)\mathit{u}=\mathit{F},\hfill \\ V^{\mathrm{lattice}}\left(\mathit{\gamma }\right)\le {\widehat{V}}^{\mathrm{lattice}},\hfill \\ l_b\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\right.\hfill \end{array}\right.$$

(11)

where $\mathit{d}$ is the static displacement vector, $\mathit{u}$ contains the magnitudes of the dynamic response, $\mathit{F}$ contains the external force amplitudes, $\mathit{K}\mathit{d}=\mathit{F}$ is the static equilibrium equation ($\Omega =0$), $(-{\Omega }^2\mathit{M}+\mathit{K})\mathit{u}=\mathit{F}$ is the dynamic equilibrium equation. $\Omega$ is the angular frequency of the time dependent harmonic load $\mathit{F}cos(\Omega t)$, and ${\widehat{V}}^{\mathrm{macro}}$ and ${\widehat{V}}^{\mathrm{lattice}}$ are the upper limits on the macro-layout volume and the lattice volume, respectively. Multiple load cases are included in the game by considering weighted compliance in (10) and weighted dynamic compliance in (11).

The players’ strategies in (10) and (11) are solved using SLP, using the sensitivities and filters presented below. The sensitivities of $c_s$ and $c_d$, using the adjoint method, are given by

$$s_e^s={\displaystyle \frac{\partial c_s}{\partial {\rho }_e}}=-{\mathit{d}}^T{\displaystyle \frac{\partial \mathit{K}}{\partial {\rho }_e}}\mathit{d},$$

(12)

$$s_e^d={\displaystyle \frac{\partial c_d}{\partial {\gamma }_e}}=-\left({\mathit{F}}^T\mathit{u}\right){\mathit{u}}^T\left(-{\Omega }^2{\displaystyle \frac{\partial \mathit{M}}{\partial {\gamma }_e}}+{\displaystyle \frac{\partial \mathit{K}}{\partial {\gamma }_e}}\right)\mathit{u},$$

(13)

where

$${\displaystyle \frac{\partial \mathit{K}}{\partial {\rho }_e}}={\displaystyle \frac{n+1}{{(1+n(1-{\rho }_e))}^2}}{\mathit{k}}_e\left({\gamma }_e\right),$$

(14)

$${\displaystyle \frac{\partial \mathit{M}}{\partial {\gamma }_e}}={\displaystyle \frac{{\rho }_e}{1+exp(-{\lambda }_m({\rho }_e-{\widehat{\rho }}_m))}}{\mathit{m}}_e$$

(15)

and

$${\displaystyle \frac{\partial \mathit{K}}{\partial {\gamma }_e}}={\displaystyle \frac{{\rho }_e}{1+n(1-{\rho }_e)}}{\displaystyle \frac{\partial {\mathit{k}}_e}{\partial {\gamma }_e}},$$

(16)

where $\partial f_{ij}/\partial {\gamma }_e$ are needed in (16). The latter is obtained by taking the derivative of (26).

The sensitivities of the volumes $V^{\mathrm{macro}}$ and $V^{\mathrm{lattice}}$ are, of course,

$${\displaystyle \frac{\partial V^{\mathrm{macro}}}{\partial {\rho }_e}}={\displaystyle \frac{\partial V^{\mathrm{lattice}}}{\partial {\gamma }_e}}=V_e.$$

(17)

The sensitivities are treated by a density filter for both ${\rho }_e$ and ${\gamma }_e$ [33], and, in addition, ${\rho }_e$ is passing a smooth Heaviside function [34]. Thus, ${\rho }_e$ is treated by

$${\rho }_e^{\mathrm{filt}}={\displaystyle \frac{{\displaystyle \sum _{g=1}^{n_{\mathrm{el}}}{\delta }_g{V}_g{\rho }_g}}{{\displaystyle \sum _{f=1}^{n_{\mathrm{el}}}{\delta }_f{V}_f}}},$$

(18)

where $n_{\mathrm{el}}$ is the number of finite elements,

$${\delta }_f={\delta }_f\left(e\right)={\left(r_{\min }-\mathrm{dist}(e,f)\right)}_{+},$$

(19)

$\mathrm{dist}(e,f)$ denotes the distance between the center of element e and f, and $r_{\min }$ is the filter radius. Equation (18) is also applied on ${\gamma }_e$. Furthermore,

$${\rho }_e^{\mathrm{heav}}={\displaystyle \frac{tanh\left(\beta \eta \right)+tanh\left(\beta ({\rho }_e^{\mathrm{filt}}-\eta )\right)}{tanh\left(\beta \eta \right)+tanh\left(\beta \right(1-\eta \left)\right)}},$$

(20)

where $\eta$ defines the threshold and $\beta$ sets the slope of the smooth Heaviside filter. In the numerical examples, $\eta =0.5$, the filter is activated after 100 SLP iterations, and then $\beta$ is ramped from 1 to 20.

4. The Multi-Scale Two-Player Game

The non-cooperative multi-scale two-player game is to find an equilibrium point $({\mathit{\rho }}^{*},{\mathit{\gamma }}^{*})$ that is optimal simultaneously for players $P_1$ and $P_2$ in (10) and (11). In this work, we try to find such equilibrium points by applying a sequential Gauss–Seidel algorithm with SLP. A recent paper on solving a two-player game in topology optimization using a Gauss–Seidel approach is presented by Thore et al. [25]. Generalized Nash equilibrium problems solved using the Gauss–Seidel method are studied in detail by Nie et al. [35].

The corresponding LP-problems for player $P_1$ and player $P_2$ in (10) and (11) read

$$P_1\left\{\begin{array}{c}\mathrm{Given}{\mathit{\gamma }}^k:\hfill \\ {\displaystyle \underset{\mathit{\rho }}{min}\sum _e{s}_e^s{\rho }_e}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{\displaystyle \sum _e{V}_e{\rho }_e\le {\widehat{V}}^{\mathrm{macro}},}\hfill \\ {\mathit{\rho }}^k-{\mathit{\delta }}^{\rho }\le \mathit{\rho }\le {\mathit{\rho }}^k+{\mathit{\delta }}^{\rho },\hfill \end{array}\right.\hfill \end{array}\right.$$

(21)

$$P_2\left\{\begin{array}{c}\mathrm{Given}{\mathit{\rho }}^k:\hfill \\ {\displaystyle \underset{\mathit{\gamma }}{min}\sum _e{s}_e^d{\gamma }_e}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{\displaystyle \sum _e{V}_e{\gamma }_e\le {\widehat{V}}^{\mathrm{lattice}},}\hfill \\ {\mathit{\gamma }}^k-{\mathit{\delta }}^{\gamma }\le \mathit{\gamma }\le {\mathit{\gamma }}^k+{\mathit{\delta }}^{\gamma },\hfill \end{array}\right.\hfill \end{array}\right.$$

(22)

where ${\mathit{\delta }}^{\rho }$ and ${\mathit{\delta }}^{\gamma }$ contain the move limits.

A strategy of conservative or aggressive move limits depending on the solution history is adopted. The strategy can be formulated as

$${\delta }_e^{\rho }=\left\{\begin{array}{cc}0.1\left(\mathrm{agressive}\right),\hfill & \mathrm{if} \Delta {\rho }_e^k\Delta {\rho }_e^{k-1}>0,\hfill \\ 0.025\left(\mathrm{conservative}\right),\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(23)

where $\Delta {\rho }_e^k={\rho }_e^{k-1}-{\rho }_e^{k-2}$. The same strategy is adopted for ${\mathit{\delta }}^{\gamma }$. In the initialization of the algorithm, $\Delta {\rho }_e^k=\Delta {\gamma }_e^k=0$ for $k<3$, implying that the three first iterations are conservative for all densities.

The sequential Gauss–Seidel algorithm reads

Set $k=0$ and guess an initial state (${\mathit{\rho }}^0,{\mathit{\gamma }}^0$), and then solve the following steps sequentially until convergence:

• Let player $P_1$ solve (21), and let ${\mathit{\rho }}^{k+1}$ be equal to the optimal solution.
• Set ${\mathit{\rho }}^k={\mathit{\rho }}^{k+1}$, let player $P_2$ solve (22), and let ${\mathit{\gamma }}^{k+1}$ be equal to the optimal solution.
• Let k = k + 1, and go to Step 1.

After a certain number of loops, the algorithm is halted and the convergence is checked. If a solution is believed to be close to an equilibrium point, then the algorithm is stopped. Otherwise, the algorithm is restarted from $({\mathit{\rho }}^k,{\mathit{\gamma }}^k)$ and continues until a new check of convergence is performed.

The proposed new non-cooperative multi-scale two-player game is implemented in our in-house toolbox using MATLAB and Fortran for three-dimensional problems with non-uniform tetrahedral meshes. This toolbox is utilized in the next section for demonstrating the numerical performance of the new game formulation.

5. Numerical Examples

The sequential Gauss–Seidel algorithm with SLP presented in the previous section is first demonstrated by solving Harker’s benchmark from [36], where the strategies of player $P_1$ and player $P_2$ read

$$P_1\left\{\begin{array}{c}\mathrm{Given}{x}_2^{*}:\hfill \\ {\displaystyle \underset{x_1}{min}{x}_1^2+{\displaystyle \frac{8}{3}}{x}_1{x}_2-34x_1}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{x}_1+x_2\le 15,\hfill \\ 0\le x_1\le 10,\hfill \end{array}\right.\hfill \end{array}\right.P_2\left\{\begin{array}{c}\mathrm{Given}{x}_1^{*}:\hfill \\ {\displaystyle \underset{x_2}{min}{x}_2^2+{\displaystyle \frac{5}{4}}{x}_1{x}_2-24.25x_2}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{x}_1+x_2\le 15,\hfill \\ 0\le x_2\le 10,\hfill \end{array}\right.\hfill \end{array}\right.$$

(24)

using variables $x_1$ and $x_2$. By performing a linearization at a state $(x_1^k,x_2^k)$, the corresponding LP-problems of (24) can be derived:

$$P_1\left\{\begin{array}{c}\mathrm{Given}{x}_2^k:\hfill \\ {\displaystyle \underset{x_1}{min}(2x_1^k+{\displaystyle \frac{8}{3}}{x}_2^k-34)x_1}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{x}_1+x_2^k\le 15,\hfill \\ 0\le x_1\le 10,\hfill \\ x_1^k-\delta \le x_1\le x_1^k+\delta ,\hfill \end{array}\right.\hfill \end{array}\right.P_2\left\{\begin{array}{c}\mathrm{Given}{x}_1^k:\hfill \\ {\displaystyle \underset{x_2}{min}(2x_2^k+{\displaystyle \frac{5}{4}}{x}_1^k-24.25)x_2}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{x}_1^k+x_2\le 15,\hfill \\ 0\le x_2\le 10,\hfill \\ x_2^k-\delta \le x_2\le x_2^k+\delta .\hfill \end{array}\right.\hfill \end{array}\right.$$

(25)

The two-player game in (24) is now solved by adopting the sequential Gauss–Seidel algorithm with SLP presented at the end of the previous section using the LP-problems in (25). Letting the initial guess be $(x_1^0,x_2^0)=(5,5)$, the convergence of the algorithm towards the equilibrium point $(x_1^{*},x_2^{*})=(5,9)$ reported by Harker [36] is presented in Figure 3.

The first benchmark demonstrating the proposed game between the players in (10) and (11) is the pin-jointed beam presented already in Figure 1. The problem has a 2D character like the following two benchmarks considered below, but these are fully treated as 3D problems using a thickness of 10 mm. Linear tetrahedral elements with a characteristic length of 2.5 mm are used and the filter radius is set to 7.5 mm. This setting is used for the three first benchmarks. In addition, for these benchmarks, the macro-volume fraction is 0.5 and the fraction of lattice in the macro-volume is 0.4, implying a total volume fraction of bulk material equals to 0.2.

Lower and upper bounds for ${\gamma }_e$ are $l_b$ = 0.2 and $u_b$ = 0.6, respectively, for all benchmarks. This setting promotes good-quality 3D printing of the designs and removal of powder to be easily done.

The number of elements of the beam is 41,753. The downward force is $F=$${10}^4$ N, and a small horizontal force of $0.1F$ is also added as a second load case, see Figure 4. The two load cases are treated using weighted compliances with weights equal to one. The extra horizontal load case is just added to prevent critical modes from being developed in this direction. For the same reason, the out-of-plane freedom for this node is constrained. Similar conditions are applied to the next two benchmarks.

Isotropic bulk properties with Young’s modulus equal to 2.1 $\times {10}^5$ MPa and Poisson’s ratio set to 0.3 are assumed for all benchmarks. The TPMS-based lattice properties are governed by the shell-based Schwarz-D structure (The choice of the Schwarz-D lattice is based on experience: it prints well, and powder is easily removed. The experience is based on the FDM, SLS, and SLM manufacturing of PLA, PA11, AlSi10Mg, 316L, and Nickel alloys. The Gyroid lattice is another good choice) using material interpolation laws $f_{ij}$ derived by Strömberg [8] (see also Figure 5). Explicitly these functions are given as the following convex combination:

$$\begin{array}{c}{f}_{ij}={\displaystyle \frac{{\alpha }_1^{ij}{\gamma }_e}{1+(1-{\gamma }_e)}}+{\displaystyle \frac{{\alpha }_2^{ij}{\gamma }_e}{1+2(1-{\gamma }_e)}}+{\displaystyle \frac{{\alpha }_3^{ij}{\gamma }_e}{1+3(1-{\gamma }_e)}}+\dots \\ {\displaystyle \frac{{\alpha }_4^{ij}{\gamma }_e}{1+5(1-{\gamma }_e)}}+{\alpha }_1^{ij}{\gamma }_e+{\alpha }_6^{ij}{\gamma }_e^2+{\alpha }_7^{ij}{\gamma }_e^3,\end{array}$$

(26)

where ${\alpha }_k^{ij}\ge 0$, ${\sum }_k{\alpha }_k^{ij}=1$ are regression coefficients presented in

Table 1. The regression coefficients for the material interpolation laws $f_{ij}$ in (26) for the shell-based Schwarz-D lattice structure.

${\mathit{f}}_{\mathit{ij}}$	${\mathit{\alpha }}_1^{\mathit{ij}}$	${\mathit{\alpha }}_2^{\mathit{ij}}$	${\mathit{\alpha }}_3^{\mathit{ij}}$	${\mathit{\alpha }}_4^{\mathit{ij}}$	${\mathit{\alpha }}_5^{\mathit{ij}}$	${\mathit{\alpha }}_6^{\mathit{ij}}$	${\mathit{\alpha }}_7^{\mathit{ij}}$
$f_{11}$	0	0	0.1837	0	0.3157	0	0.5006
$f_{21}$	0	0	0.7155	0	0.0762	0	0.2083
$f_{44}$	0	0	0	0	0.2429	0.5650	0.1921

.

The optimal grading plotted on the macro-layout of the beam for the pure static load case ($\Omega =0$) is presented in Figure 4 (the static design). The corresponding stl-file is presented in Figure 1. Modal analysis with structural damping of this design is performed using Abaqus/Standard and mode 5 is identified as critical with an eigenfrequency at 1604 Hz. The first four modes are out-of-plane modes. The structural damping factor is set to 0.02 for all benchmarks in order to bind the frequency responses. A reminder on structural damping is presented in Appendix A.

Next, the pin-jointed beam is optimized by letting the excitation frequency $\Omega$ be equal to 1604 Hz (The same symbol $\Omega$ is used for the harmonic angular frequency (rad/s) and the harmonic frequency (Hz)). Resonance will of course occur for the static design for this excitation. The optimal grading on the macro-layout for the dynamic design obtained by solving the game defined by (10) and (11) for $\Omega =1604$ Hz is presented in Figure 4. The corresponding stl-file is also presented in the same figure. The macro-layout is similar to the Static one, but the grading differs significantly. By zooming in on the figure, the grading of the lattice also appears clearly for the stl-file in gray. Theoretically, the lattice scale shown in Figure 4 should be significantly smaller. However, the scale used here, and in the subsequent analyses, is selected with practical considerations for 3D printing in mind. It remains sufficiently small to accurately capture the lowest eigenmodes as well as the global stiffness.

The new frequency of the critical mode 5 is shifted from 1604 to 1178 Hz. Furthermore, a large dashed window of frequencies where the performance of the dynamic design is better than the static one is also identified. Here, and in the following, the window of better performance means that the displacement amplitude obtained for the dynamic design is smaller than the corresponding amplitude for the static design.

The convergence of the static and the dynamic compliances during the game is presented in Figure 3. It is clear that the game converges to an equilibrium point. This behavior is observed for all benchmarks. The activation of the Heaviside filter presented in (20) after 100 iterations also appears clearly in this convergence plot. Thus, the computational performance is good, and the total computational cost is only slightly larger than twice the cost of the corresponding single static compliance problem.

The second benchmark is a clamped beam with boundary conditions according to Figure 6. The design domain is meshed using 41,350 linear tetrahedral elements. The optimal macro-layout with the grading of the static design ($\Omega =0$) is presented in the figure. The eigenfrequency of the critical mode 3 of this static design is 2382 Hz. A dynamic design for $\Omega =2383$ Hz is then generated using the game formulation and it is also presented in the figure. The eigenfrequency of the critical mode 3 is now shifted from 2383 Hz to 1799 Hz. A large dashed window with frequencies where the new dynamic design performs better than the Static one is also presented. The macro-layouts are slightly different, but the grading differs significantly between the static and dynamic design. An stl-file of the dynamic design is also depicted, where the grading appears clearly by zooming.

The standard L-shaped benchmark is considered in Figure 7. A total of 50,486 linear tetrahedral elements are used for the design domain. The macro-layout with grading for the static design is presented in the figure. In Figure 8 the modal analysis for this design is presented. Two critical modes are identified, Mode 2 at the lower frequency of 372 Hz and Mode 6 at the higher frequency of 1206 Hz. Two corresponding optimal dynamic designs are then generated using the game formulation: Dynamic-L for $\Omega =372$ Hz and Dynamic-H for $\Omega =1206$ Hz. These designs are presented in Figure 7, both as density gradings on the macro-layouts as well as corresponding stl-files.

The modal analyses of the optimal designs are presented in Figure 8. The lower critical frequency is shifted from 372 Hz to 443 Hz (upwards), and the higher critical frequency is shifted from 1206 Hz to 859 Hz. Dashed windows where the dynamic designs perform better than the static design are also identified. A new cross member is developed for the macro-layout of the Dynamic-L design and the grading turns to the lower limit at the applied force. For the Dynamic-H design, the macro-layout is quite similar to the static design, but the grading differs significantly.

Finally, a three-dimensional L-shaped benchmark similar to the one presented by Borrvall and Petersson [31] is considered in Figure 9. This benchmark was also considered in [7]. The design domain is discretized using 223,471 linear tetrahedral elements. The filter radius is set to 10 mm. The macro-volume fraction is 0.25 and the fraction of lattice in the macro-volume is 0.4, implying a total volume fraction of bulk material equals to only 0.1 (The standard compliance problem with this low volume fraction of 0.1, and identical settings in mesh, filter radius and load cases does not converge to a proper solid design). The four corners of the back-side are fixed, and the lower opposite corner is subjected to a weighted load case using three orthogonal forces F with equal magnitudes and equal weights.

First, an optimal static design is obtained by letting $\Omega$ = 0. This design is plotted for the first fundamental mode with an eigenfrequency of 986 Hz. Here, 2,045,823 elements linear tetrahedral elements are used in the modal analysis. Then, a dynamic design is generated by the game for $\Omega$ = 986 Hz. The optimal dynamic design is presented for Mode 2 with an eigenfrequency of 759 Hz. This is the shift of the second eigenfrequency of the static design at 1078 Hz. The first eigenfrequency of the static design at 986 Hz is shifted to 678 Hz. A well-defined dashed window with frequencies at the critical frequency of 986 Hz is obtained where the dynamic design is better than the Static one. In addition, a dashed window for higher frequencies with a similar trend also appears due to the fact that the critical frequency at 1500 Hz (Mode 3) for the static design is shifted to 1164 Hz (Mode 3) for the new dynamic design. In order to enlarge the size of the dashed windows, one might suggest using multiple load cases with multiple harmonic frequencies. Then, each load case $\mathit{F}$ is considered for N frequencies ${\Omega }_i$ ($n=1,\dots ,N$) such that (11) modifies to

$$P_2^M\left\{\begin{array}{c}\mathrm{Given}{\mathit{\rho }}^{*}:\hfill \\ {\displaystyle \underset{(\mathit{\gamma },\mathit{u})}{min}{c}_d=\sum _{i=1}^N{\displaystyle \frac{1}{2}}{\left({\mathit{F}}^T{\mathit{u}}_i\right)}^2}\hfill \\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\left(-{\Omega }_i^2\mathit{M}({\mathit{\rho }}^{*},\mathit{\gamma })+\mathit{K}({\mathit{\rho }}^{*},\mathit{\gamma })\right){\mathit{u}}_i=\mathit{F},(i=1,\dots ,N)\hfill \\ V^{\mathrm{lattice}}\left(\mathit{\gamma }\right)\le {\widehat{V}}^{\mathrm{lattice}},\hfill \\ l_b\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1}.\hfill \end{array}\right.\hfill \end{array}\right.$$

(27)

Considering multiple harmonic frequencies is a topic for future work.

6. Concluding Remarks

In this paper, we propose, implement, and demonstrate a new novel multi-scale non-cooperative two-player game for avoiding resonance under specified harmonic loads. Player 1 minimizes static compliance, while Player 2 minimizes dynamic compliance. The game is solved using a Gauss–Seidel algorithm with sequential linear programming. It is numerically demonstrated that equilibrium designs with optimal macro-layouts featuring optimal local grading of shell-based Schwarz-D lattice structures are generated with natural eigenfrequencies significantly different from the excitation frequency of the harmonic loads. Additionally, the proposed multi-scale game effectively prevents the development of designs with rigid body modes, a common issue when static compliance is not constrained in standard dynamic compliance formulations. This constraint is managed by Player 1 in the proposed game. In the near future, it would be intriguing to adopt this multi-scale game approach for other types of topology optimization problems.