Mixed-integer second-order cone optimization for composite discrete ply-angle and thickness topology optimization problems

Abstract

Discrete variable topology optimization problems are usually solved by using solid isotropic material with penalization (SIMP), genetic algorithms (GA), or mixed-integer nonlinear optimization (MINLO). In this paper, we propose formulating discrete ply-angle and thickness topology optimization problems as a mixed-integer second-order cone optimization (MISOCO) problem. Unlike SIMP and GA methods, MISOCO efficiently finds the problem’s globally optimal solution. Furthermore, in contrast with existing MISOCO formulations of discrete ply-angle optimization problems, our reformulations allow the structure to change topology, consider the more realistic Tsai–Wu stress yield criteria constraint, and eliminate checkerboard patterns using simple linear constraints. We address two types of discrete ply-angle and thickness problems: a structural mass minimization problem and a compliance optimization problem where the objective is to maximize the structural stiffness. For each element, one first chooses if the element is present or not in the structure. One can then choose the element’s ply-angle and thickness from a finite set of possibilities for the former case. The discrete design space for ply-angle and thickness is a result of manufacturing limitations. To improve the problem’s MISOCO solution approach, we develop valid inequality constraints to tighten the continuous relaxation of the MISOCO reformulation. We compare the performance of various MISOCO solvers: Gurobi, CPLEX, and MOSEK to solve the MISOCO reformulation. We also use BARON to solve the original MINLO formulations of the problems. Our results show that solving the MISOCO problem’s formulation using MOSEK is the most efficient solution approach.

Appendix

1.1 Constitutive matrix calculation

1.1.1 In-plane constitutive relation

In this section, we present the constitutive matrix ${\mathbf {Q}}_i$ derivation procedures. In a frame that the x axis is aligned with the fiber direction, we have the following constitutive relationships:

$$\begin{aligned} \begin{bmatrix} \sigma _1\\ \sigma _2\\ \tau _{12} \end{bmatrix} = {\mathbf {Q}} \begin{bmatrix} \epsilon _1\\ \epsilon _2\\ \epsilon _{12} \end{bmatrix}, \end{aligned}$$

where ${\mathbf {Q}}$ is given as

$$\begin{aligned} {\mathbf {Q}}= \begin{bmatrix} Q_{11} &{} \quad Q_{12} &{} \quad 0\\ Q_{21} &{} \quad Q_{22} &{} \quad 0\\ 0 &{} \quad 0 &{} \quad 2Q_{66}\\ \end{bmatrix}, \end{aligned}$$

and

$$\begin{aligned} Q_{11}& = \frac{E_1}{1-\nu _{21}\nu _{12}}\\ Q_{12}& = \frac{\nu _{12} E_2}{1-\nu _{21}\nu _{12}}=Q_{21}\\ Q_{22}& = \frac{E_2}{1-\nu _{21}\nu _{12}}\\ Q_{66}& = G_{12}. \end{aligned}$$

For more general cases where the x axis and the fiber direction are with an angle of $\theta$, using the laws of tensor transformation relations, we have:

$$\begin{aligned} \begin{bmatrix} \sigma _x\\ \sigma _y\\ \tau _{xy} \end{bmatrix} = {\mathbf {T}}^{-1}(\theta ) \begin{bmatrix} \sigma _1\\ \sigma _2\\ \tau _{12} \end{bmatrix}, \end{aligned}$$

where

$$\begin{aligned} {\mathbf {T}}^{-1}(\theta )= \begin{bmatrix} \cos ^2\theta &{} \quad \sin ^2\theta &{} \quad -2\cos \theta \sin \theta \\ \sin ^2\theta &{} \quad \cos ^2\theta &{} \quad 2\cos \theta \sin \theta \\ \cos \theta \sin \theta &{} \quad -\cos \theta \sin \theta &{} \quad \cos ^2\theta - \sin ^2\theta \end{bmatrix}, \end{aligned}$$

and

$$\begin{aligned} \begin{bmatrix} \epsilon _1\\ \epsilon _2\\ \epsilon _{12} \end{bmatrix} = {\mathbf {T}}(\theta ) \begin{bmatrix} \epsilon _x\\ \epsilon _y\\ \epsilon _{xy} \end{bmatrix}. \end{aligned}$$

Thus, the general constitutive law is given as:

$$\begin{aligned} \begin{bmatrix} \sigma _x\\ \sigma _y\\ \tau _{xy} \end{bmatrix} = {\mathbf {T}}^{-1}(\theta ) {\mathbf {Q}} {\mathbf {T}}(\theta ) \begin{bmatrix} \epsilon _x\\ \epsilon _y\\ \epsilon _{xy} \end{bmatrix}. \end{aligned}$$

For the jth choice with the angle as $\theta _j$, we define:

$$\begin{aligned} {\mathbf {Q}}_j={\mathbf {T}}^{-1}(\theta _j) {\mathbf {Q}} {\mathbf {T}}(\theta _j). \end{aligned}$$

We can give an explicit expression for ${\mathbf {Q}}_j$ by expanding the previous equation:

$$\begin{aligned} Q_{j, 11}& = Q_{11}\cos ^4\theta _j + Q_{22}\sin ^4\theta _j+2(Q_{12}+2Q_{66})\sin ^2\theta _j\cos ^2\theta _j\\ Q_{j, 12}& = (Q_{11}+Q_{22}-4Q_{66})\cos ^2\theta _j\sin ^2\theta _j+Q_{12}(\sin ^4\theta _j+\cos ^4\theta _j)\\ Q_{j, 22}& = Q_{11}\sin ^4\theta _j + Q_{22}\cos ^4\theta _j+2(Q_{12}+2Q_{66})\sin ^2\theta _j\cos ^2\theta _j\\ Q_{j, 16}& = (Q_{11}-Q_{12}-2Q_{66})\cos ^3\theta _j\sin \theta _j-(Q_{22}-Q_{12}-2Q_{66})\cos \theta _j\sin ^3\theta _j\\ Q_{j, 26}& = (Q_{11}-Q_{12}-2Q_{66})\cos \theta _j\sin ^3\theta _j-(Q_{22}-Q_{12}-2Q_{66})\cos ^3\theta _j\sin \theta _j\\ Q_{j, 66}& = (Q_{11}+Q_{22}-2Q_{12}-2Q_{66})\cos ^2\theta _j\sin ^2\theta _j+Q_{66}(\sin ^4\theta _j+\cos ^4\theta _j).\\ \end{aligned}$$

1.1.2 Out-of-plane constitutive relation

The out-of-plane constitutive relationship is given as

$$\begin{aligned} \begin{bmatrix} \sigma _{yz}\\ \sigma _{xz} \end{bmatrix} = {\mathbf {C}}_j \begin{bmatrix} \psi _y + \frac{\partial w}{\partial y}\\ \psi _x + \frac{\partial w}{\partial x}, \end{bmatrix}, \end{aligned}$$

where

$$\begin{aligned} {\mathbf {C}}_j= \begin{bmatrix} C_{j, 44} &{} \quad C_{j, 45}\\ C_{j, 45} &{} \quad C_{j, 55} \end{bmatrix}, \end{aligned}$$

and

$$\begin{aligned} C_{j, 44}& = C_{44}\cos ^2\theta _j + C_{55}\sin ^2\theta _j\\ C_{j, 45}& = (C_{55} - C_{44})\cos \theta _j\sin \theta _j, \\ C_{j, 55}& = C_{44}\sin ^2\theta _j + C_{55}\cos ^2\theta _j. \end{aligned}$$

We also know that $C_{44}=G_{23}$ and $C_{55}=G_{12}$.

1.2 Strain operator

We calculate the strain of the upper surface center for each element that is later used in Tsai–Wu failure criterion. The in-plane strain is given as:

$$\begin{aligned} \begin{aligned} \pmb {\epsilon }&= \pmb {\epsilon }_l + z{\chi }, \\ \pmb {\epsilon }_l&={\mathbf {L}}_1{\mathbf {u}}, \\ {\chi }&={\mathbf {L}}_2{\mathbf {u}}, \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} {\mathbf {L}}_1&=\left[ \begin{array}{c|ccccc|c} &{} \quad N_{i,x} &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad \\ \dots &{} \quad 0 &{} \quad N_{i,y} &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad \dots \\ &{} \quad N_{i,y} &{} \quad N_{i,x} &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad \\ \end{array}\right] ,\\ {\mathbf {L}}_2&=\left[ \begin{array}{c|ccccc|c} &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad N_{i,x} &{} \quad 0 &{} \quad \\ \dots &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad N_{i,y} &{} \quad \dots \\ &{} \quad 0 &{} \quad 0 &{} \quad 0 &{} \quad N_{i,y} &{} \quad N_{i,x} &{} \quad \\ \end{array}\right] .\\ \end{aligned} \end{aligned}$$

For the shape function $N_i$, we have

$$\begin{aligned} \begin{aligned} N_{i+3j-3}&= f_i(\eta ) f_j(\xi ), \,\, i,j=1,2,3,\\ f_1(x)&= \frac{1}{2}x(x-1),\\ f_2(x)&= 1-x^2,\\ f_3(x)&= \frac{1}{2}(1+x)x. \end{aligned} \end{aligned}$$

To evaluate the strain of the center point over the top surface, we simply set $(\eta , \xi , z)=(0, 0, h/2)$. Finally,

$$\begin{aligned} \pmb {\epsilon }_i={\mathbf {B}}({\mathbf {x}}_i){\mathbf {u}}_i =\left( {\mathbf {L}}_1(0)+\frac{h}{2}{\mathbf {L}}_2(0)\right) {\mathbf {u}}_i. \end{aligned}$$

1.3 Local stiffness matrix calculation

In this section, we present the derivation of the local stiffness matrix ${\mathbf {K}}_{ij}$ derivation, which is based on classic lamination theory (CLT). We begin the derivation by giving the expression for matrices ${\mathbf {A}}, {\mathbf {B}}, {\mathbf {D}}$ representing the extensional stiffness matrix, the coupling stiffness matrix and the bending stiffness matrix, respectively:

$$\begin{aligned} \begin{bmatrix} {\mathbf {A}}&\quad {\mathbf {B}}&\quad {\mathbf {D}} \end{bmatrix} = \int _{-h_i/2}^{h_i/2} \begin{bmatrix} 1&\quad z&\quad z^2 \end{bmatrix} {\mathbf {Q}}_j dz. \end{aligned}$$

The traverse shear stiffness matrix is given as:

$$\begin{aligned} {\mathbf {S}}=\kappa \int _{-h_i/2}^{h_i/2}{\mathbf {C}}_j dz, \end{aligned}$$

where $\kappa =\frac{5}{6}$. Notice that to be consistent with the assumption made in this paper here we assume that there is only one layer of material.

We now present the stiffness matrix expression:

$$\begin{aligned} \begin{aligned} {\mathbf {K}}_{ij}&= {\mathbf {K}}_{ij,L_1} + {\mathbf {K}}_{ij,L_2}\\&=\int _{A_i} \left( L_1^\intercal {\mathbf {A}} L_1 +L_1^\intercal {\mathbf {B}}L_2 + L_2^\intercal {\mathbf {B}}L_1 + L_2^\intercal {\mathbf {D}} L_2 \right) dA_i + \int _{A_i} \left( L_3^\intercal {\mathbf {S}}L_3\right) dA_i \end{aligned} \end{aligned}$$

Then we apply the Gaussian quadrature for numerical integration. For $K_{ij,L_1}$, we apply $3\times 3$ integration points; while for $K_{ij,L_2}$, we apply $2\times 2$ integration points. This helps to avoid “shear locking”.

1.4 A brief overview of SOC

A SOC constraint can be expressed in the following format (Lobo et al. 1998)

$$\begin{aligned} \Vert {\mathbf {A}} {\mathbf {x}} + {\mathbf {b}} \Vert \le {\mathbf {c}}^\top {\mathbf {x}} + d, \end{aligned}$$

(27)

where $\Vert {\cdot }\Vert$ denotes the Euclidean norm, ${\mathbf {x}}\in {{\mathbb {R}}}^n$ is the variable vector, ${\mathbf {A}} \in {{\mathbb {R}}}^{m \times n}$, ${\mathbf {b}}\in {{\mathbb {R}}}^{m}$, ${\mathbf {c}}\in {{\mathbb {R}}}^n$, and $d\in {{\mathbb {R}}}$.

A SOC constraint appears naturally when the feasible set of a convex quadratic function is considered. In particular, let ${\mathbf {Q}}$ be PSD, then

$$\begin{aligned} \{(y, {\mathbf {x}})\in {{\mathbb {R}}}\times {{\mathbb {R}}}^n: y^2 \ge {\mathbf {x}}^\intercal {\mathbf {Q}} {\mathbf {x}}\}, \end{aligned}$$

(28)

is equal to the following set

$$\begin{aligned} \{(y, {\mathbf {x}})\in {{\mathbb {R}}}\times {{\mathbb {R}}}^n: \Vert {\mathbf {Q}}^{1/2}{\mathbf {x}}\Vert \le y\}. \end{aligned}$$

(29)

Many more examples, including truss topology optimization problems, showing how SOC constraints arise in practice are discussed by Lobo et al. (1998).