## 2 Design of Base Cells
### Topology Optimization and Sensitivity Analysis  
Combining the PLSM and numerical homogenization theory, the topology optimization model for cellular structure with extreme properties can be established:  
\[find: \alpha_{i}(i=1,2,3...N) \tag{7}\] \[Maximize:J(u,\Phi)=\sum\limits_{i,j,k=1}^{d}\alpha_{ijkl}D_{ijkl }^{H}(u,\Phi)\] \[Subject\ to:\left\{ \begin{array}{c}G(\Phi)=\int_{V}H(\Phi)dV \leq V_{\max}\\ a(u,v,\Phi)=l(u,v)\\ \alpha_{\min}\leq\alpha_{i}\leq\alpha_{\max}\end{array} \right.\]  
In the optimization model, \(\alpha_{i}\) is the expansion coefficient of level set function \(\Phi\), which works as the design variable. \(\alpha_{\min}\) and \(\alpha_{\max}\) respectively represent the lower and upper bounds of the design variables. In the PLSM, \(\alpha_{\min}\) can be selected as \(2\times\min(\alpha_{i}^{(\gamma-1)})\), \(\alpha_{\max}\) can be defined as \(2\times\max(\alpha_{i}^{(\gamma-1)})\)[51]. Here, \(\gamma\) represent the iteration number, and \(\alpha_{i}^{(\gamma-1)}\) denotes the design variables in the last iteration. \(N\) is the total number of design variables, which is also the number of the level set nodes or finite element nodes [50; 51]. \(\omega_{ijkl}\) in the objective function represents the weight factors of different parameters in the elasticity tensor \(D^{H}\). \(G(\Phi)\) is the volume constraint with a maximum volume \(V_{\max}\). \(H(\Phi)\) is the Heaviside function used to distinguish different parts in the design domain. Considering a smooth Heaviside function, \(H(\Phi)\) can be defined as [20]:  
\[H(\Phi)= \begin{cases}\eta,&\Phi<-\Delta\\ \frac{3(1-\eta)}{4}(\frac{\Phi}{\Delta}-\frac{\Phi^{3}}{3\Delta^{3}})+\frac{1+ \eta}{2},&-\Delta\leq\Phi\leq\Delta\\ 1,&\Phi>\Delta\end{cases}  \tag{8}\]  
In case of the numerical singularity, \(\eta\) is set to be a small positive number. \(\Delta\) represents the width of the numerical approximation. The equation of bilinear energy form and linear load form are respectively defined as follows:  
\[\left\{ \begin{array}{c}a(u,v,\Phi)=\int_{V}\varepsilon_{pq}^{*}(u^{(ij)})D_ {pqrs}\varepsilon_{rs}^{*}(v^{(ij)})H(\Phi)dV\\ l(v,\Phi)=\int_{V}\varepsilon_{pq}^{0(ij)}D_{pqrs}\varepsilon_{rs}^{*}(v^{(ij)}) H(\Phi)dV\end{array} \right. \tag{9}\]To solve the optimization Formulation (7), the derivatives of the objective function and the volume constraint with respect to the design variables should be computed. Based on the theory of shape derivative [21], the sensitivity analysis can be conducted, and the derivatives of \(D^{H}(u,\Phi)\)\(\mu(u,v,\Phi)\)\(\jmath(v,\Phi)\) with respect to time \(t\) are respectively described by:  
\[ \begin{array}{c}\frac{\partial D_{ijkl}^{H}}{\partial t}=\frac{1}{|V|}\int_{ V}(\varepsilon_{pq}^{0(ij)}-\varepsilon_{pq}^{*}(u^{(ij)}))D_{pqrs}( \varepsilon_{rs}^{0(kl)}-\varepsilon_{rs}^{*}(u^{(kl)}))V_{n}|\
abla\Phi|\delta( \Phi)dV\\ \qquad\qquad-\frac{2}{|V|}\int_{V}(\varepsilon_{pq}^{*}(\dot{u}^{(ij)}))D_{ pqrs}\Big{(}\varepsilon_{rs}^{0(kl)}-\varepsilon_{rs}^{*}(u^{(kl)})\Big{)}H(\Phi)dV \end{array}  \tag{10}\]  
\[ \begin{array}{c}\frac{\partial\dot{a}}{\partial t}=\int_{V} \varepsilon_{pq}^{*}(u^{(ij)})D_{pqrs}\varepsilon_{rs}^{*}(v^{(ij)})V_{n}| \
abla\Phi|\delta(\Phi)dV\\ \qquad\qquad+\int_{V}(\varepsilon_{pq}^{*}(\dot{u}^{(ij)})D_{pqrs}\varepsilon_ {rs}^{*}(v^{(ij)})+\varepsilon_{pq}^{*}(u^{(ij)})D_{pqrs}\varepsilon_{rs}^{*}( \dot{v}^{(ij)}))H(\Phi)dV\end{array}  \tag{11}\]  
\[ \begin{array}{c}\frac{\partial l}{\partial t}=\int_{V}\varepsilon_{pq}^{(ij) }D_{pqrs}\varepsilon_{rs}^{*}(v^{(ij)})V_{n}|\
abla\Phi|\delta(\Phi)dV+\int_{V }\varepsilon_{pq}^{(ij)}D_{pqrs}\varepsilon_{rs}^{*}(\dot{v}^{(ij)})H(\Phi)dV \end{array}