Received: 14 Jul 2018

Revised: 10 Oct 2018

Accepted: 10 Oct 2018

Published online: 11 Oct 2018

Jihong Zhu^{ 1, 2, 3, Email }Yubo Zhao^{ 1} Weihong Zhang^{ 1, 2, Email} Xiaojun Gu^{ 3} Tong Gao^{ 1, 3}Jie Kong^{ 4}Guanghui Shi^{ 5}Yingjie Xu^{ 1, 2}Dongliang Quan^{ 5}

^{1 }State IJR Center of Aerospace Design and Additive Manufacturing, School of Mechanical Engineering, Northwestern Polytechnical University, Xian, Shaanxi 710072, China

^{2 }MIIT Lab of Metal Additive Manufacturing and Innovative Design, NPU-QMUL Joint Research Institute, Northwestern Polytechnical University, Xian, Shaanxi 710072, China

^{3 }Institute of Intelligence Material and Structure, Unmanned System Technologies, Northwestern Polytechnical University, Xian, Shaanxi 710072, China

^{4 }Shaanxi Key Laboratory of Macromolecular Science and Technology, School of Science, Northwestern Polytechnical University, Xian, Shaanxi 710072, China

^{5 }Beijing Aerospace Technology Institute, Beijing 100074, China

*E-mail: JH.Zhu@nwpu.edu.cn; zhangwh@nwpu.edu.cn

A fish operculum bone is a high-efficiency load-carrying structure, which transfers the surface pressure of the fluid through the special supporting structures on the thin-walled panel to the joint position connecting the body. In order to bring this concept of bionic structures into the flight vehicle rudder structures’ optimization design, problems such as the structural feature and parameterized modeling, load-carrying behaviors analysis, optimization design method, performance evaluation etc. need to be properly solved. To this end, we studied the operculum structure and found that the properly distributed Y-shape structural branches undertake the surface load perfectly. Meanwhile, the influence of these Yshape branches distribution was discussed with numerical simulation and topology optimization. By considering the Y-shape branches as a special structural features, the bio-inspired design procedure was then established by integrating the structural layout and sizing parameters as well as their simultaneous feature-driven optimization. We also designed a typical flight vehicle rudder structure and fabricated it with stereolithography based resin additive manufacturing. With more than 20 % improvement in stiffness and strength compared with the traditional design, the advantages of the bio-inspired optimization have been clearly demonstrated.

Rudder is a key component to realize lift and maneuverability for flight vehicles. Its lightweight and high bearing capacity design is of great importance for the performance of the whole flight vehicle. Traditionally, the rudder is a typical stiffeners reinforced thin-walled panel. These stiffeners are distributed uniformly to undertake the surface load. For dozens of years, this kind of structural configuration has been widely used in flight vehicle and aerospace structure design due to its manufacturability and balance in load-carrying performances.1

In the least-weight and performance design of flight vehicle structures, sizing and shape optimizations are two traditional techniques and have been widely employed.2,3 More optimization studies were also carried out by redistributing and resizing the structural members etc. using multi-level methods.4,5 Wang et al.6 also presented a multi-stage and simultaneous layout-sizing optimization scheme for these stiffener reinforced panels, which is a significant improvement of the traditional structure pattern. Recently, due to the rapid development and severe requirements of flight vehicle industries, the weight saving and performance improvement of these regularly designed structures have almost reached the limit.7,8

In 1980s, topology optimization method was proposed9,10 and has revolutionized the traditional design patterns. By designing the material layout in the structural space, the optimal load-carrying paths from loads to boundary conditions are realized.11,12 Many successful applications have been found in flight vehicle and aerospace structures8. In these works, topology optimization achieved various sorts of beautiful and reasonable truss-like features, which have been well accepted and can be designed undertake multidisciplinary loading conditions such as static, dynamic, thermal, fluid, electromagnetic etc. with significant weight saving compared with the traditional design.

For many years, researchers have attempted to use topology optimization to design the stiffeners layout of thin-walled panels,9,13,14,15,16 especially for the flight vehicle wings, rudders, stabilizers and control surfaces. However, compared with the beautiful truss-like structures, the optimized structures here presented some patterns that are very different from the regular configurations and are difficult to be realized in engineering. Researchers also presented some topologically optimized results based on traditional flight vehicle structural configurations,17,18 which turns out to be a compromise and will lead to unsatisfied performance optimization. Alternatively, other researchers tried to optimize the stiffeners using direct layout optimization. Stiffeners reinforced panels were considered as composites where the mechanical properties of different unit cells were calculated using homogenization.19,20 Some recent developments of the stiffener designs have successfully realized curved stiffener designs.21,22

In view of this, finding proper structures from the nature as the references of structural feature and parameterized modeling for optimized results and establishing a new bio-inspired topology optimization method will be a reasonable solution. Consider the loading characteristics of the flight vehicle rudder, we chose the fish operculum bone in this paper, which is composed of some particular highefficiency load-carrying structural branches and transfers the surface loads to the joint position on body. The load-carrying behaviors of the operculum will be analyzed and the design methodology will be established accordingly

**Fig 1.** Operculum bone structures and its Y-shape branches.

**Fig 2.** Finite element model of the operculum structure and its Y-shape branches as the basic structural features.

The operculum is a facial support structure and a protective covering for the gills of a bony fish.23,24 It is also used for its breathing and feeding. Figure 1 shows the detailed operculum structure configuration of an aristichthys nobilis (bighead carp). The main load-carrying structure appears as a tree-like pattern which consists of many Y-shape branches with multi-fractal features. It is also noticed that the Y-shape branches distributed in three dimensions, which is also indicated in Fig. 1. Similar structures can be found in some concentrated force diffusion structures25,26 and thermal conduction devices.27,28

Considering these Y-shape branches as the basic structural feature of the global load-carrying structure, we chose their distribution in the panel surface, lengths on the span direction and sizes of the crosssections etc. as the parameters for numerical modelling, which can be seen in Fig. 2. The structure modelled as an assemblage of the Y-shape branches of different parameters in accordance with the real operculum structure is discretized into refined finite elements and applied with a fluid pressure load on the panel surface. The distribution of the stress level and the principal stress of one of the Y-shape branches are also presented in Fig. 3. With this load condition, compressions are mostly seen in these Y-shape branches. It is also found that the alignment of the Y-shape branches appears fine conjunctions with the directions of the principal stresses, where the surface loads converge from the panel to the branches and the root of the structure in turn, just like river flows.

To further verify the reasonableness of the structural distribution, a standard compliance minimization and density based topology optimization with the same loading condition and 10 % volume fraction for this operculum structure is carried out with the optimized material distribution presented in Fig. 4. Compared with the optimized structure where the detailed branches are not clear enough as expected, the general material distribution and the gradient panel thickness of the operculum structure are mostly alike, which indicates that the natural structure of the operculum is reasonable and optimal in surface load carrying. Moreover, the thickness of the branches as well as the thickness distribution of the panel well corresponds to the topologically optimized results.

Inspired by the above simulation and analysis, we can combine the advantages of the topologically optimized structure and the operculum

**Fig 3.** The distribution of the stress level and the principal stress of the operculum structure under a fluid pressure load on the panel surface.

**Fig 4.** Topologically optimized result of the operculum structure, where clear material distribution and gradient panel thickness of the operculum structure are obtained.

**Fig 5.** Skeleton and the modelling of the Y-shape branches of different levels, where we have branches locating their root on other branches.

structure to implement the engineering structures design. In this case, the topologically optimized structure without detailed structural branches can be employed as a reference for the initial configuration design. Meanwhile, the operculum structure consists of Y-shape branches with different patterns can be considered as a basic feature of the structure in optimization. In this way, the optimization is considered as bio-inspired feature-driven topology optimization scheme that can be formulated as follows.

(1) Design variables:

The design variables include the shape and layout parameters of the Yshape branches. As illustrated in Fig. 2, parameters in design variable vector Y are used for a single Y-shape branch to describe its geometry, location and orientation. There are 4 points, i.e. one root, two tips, one branching position and their locations defining the skeleton of the branch. At these positions, the local heights and widths are assigned to define the size of the branch. We can also have branches locating their root on other branches, which creates Y-shape branches of different levels as shown in Fig. 5. In this way, the structural branches are actually defined by a skeleton and the widths and heights at different root, branching and tip positions. The skeleton can be immediately divided into numbers of straight sections by the coordinates of the key points at the root, branching and tip positions. Each section will be modelled as a quadrangular frustum pyramid as show in Fig. 5. These shapes will be later merged together using Boolean operation. As a result, the coordinates of the key points on the skeleton, the width and height at each position will be selected as the design variables composing the vector Y.

Meanwhile, the thicknesses at different positions on the panel, which can be described as the shell element thicknesses (design variable vector T) on the panel are also introduced as the design variables. t e indicates the thickness of the eth element on the panel. (2) Objective function:

The global strain energy, which is also mentioned as the mean compliance is chosen as the design objective to be minimized. It is a compliance based topology optimization enhancing the structural global stiffness. U is the displacement vector and K is the structural stiffness matrix.

(3) Design constraints:

The first line indicates the equilibrium equation of the static finite element model. The second line defines the total volume constraint to limit the structural weight. For practical flight vehicle structure design, more constraints on the structural stress and deformations in particular positions shall be taken into account, which can be indicated as:

Here indicates the eth stress value in the structure with the total e σ number n . U is the displacement on the dth DOF with the total number e d n . The maximum displacements of a rudder always appear at the d corners. The displacement constraints are therefore defined at these positions. Then large numbers of stress constraints will lead to costly computing time when the optimization problem is solved iteratively. It is therefore crucial to use some constraint aggregation approaches such as KS function etc. The standard form of KS function is

where N stands for the number of constraints. μ is the aggregation g parameter which is assigned as 15 in this paper. Generally, the constraints should be normalized into the same scale, which will be much helpful for the constraint aggregation, especially when the constraints are in different magnitudes. More details about KS function can be seen in Martins & Poon,29 Poon & Martins30 and Gao et al.31

As the proposed bio-inspired topology optimization scheme has been defined as a simultaneous layout and thickness optimization, the definition of the design variables, i.e. the pseudo-densities in the density based topology optimization will not be applicable here. Moreover, commonly used body-fit finite element mesh will certainly lead to element remeshing in layout optimization. Consequently, we proposed to use a background mesh and a feature-driven topology optimization to undertake this problem32,33. Without loss of generality, a 2D case is illustrated as shown in Fig. 6. The original design domain will be discretized into fine and regular finite elements. Then the Y-shape branches are introduced into the design domain as numbers of features. The elements totally included inside the features are assigned as the solid material, those totally outside the features are assigned as the void material. Difficulties are found for the definition of the materials for the intermediate elements cut by the edges of the features. Three different ways are usually used to tackle the problem, i.e. the extended finite element method (XFEM34), the finite cell method (FCM35,36) and the stiffness interpolation method.14,37

XFEM locally enriched the intermediate finite element approximation by considering special treatments in shape sensitivity analysis. By defining the structural features using the level-set functions and Heaviside functions, the elastic matrices and displacement fields of these enriched elements are then precisely expressed. The FCM can be considered as a particular XFEM. It uses the high-order Ansatz functions such as Legendre polynomial to approximate the extended solution. The Gauss points of each intermediate element are locally refined in a quadtree subdivision (octree subdivision in 3D cases). XFEM and FCM hold the advantage of precise calculation for local behaviors such as displacements and stresses. However, the modelling and calculation procedure is rather complicated compared with the original finite element analysis. The stiffness interpolation method is easily implemented by interpolation the elements cut by the feature contour with intermediate stiffness between the solid and void. The stiffness value is decided by the ratio of the solid part and void part in the corresponding element. Obviously, stiffness interpolation can only concern some structural global behaviors and will sacrifice the detailed behaviors in local stresses and displacements. In this paper, as the stresses and some local displacements are taken into account in the optimization, we chose using FCM with octree subdivision.

Consider the Y-shape branches in FCM, the key to the refinement of the Gauss points is to determine the elements cut by the contours of the branches. The geometry of the branches that are modelled as quadrangular frustum pyramid will be defined as a level-set function .

**Fig 6.** Definitions of the intermediate elements cut by the edges of the features: XFEM, FCM and stiffness interpolation

Suppose there is a global structural region Ω, which is composed of

Ω is the solid region inside the structural branches with its boundary Y Г( =0). Y f A Heaviside function H is also introduced. Then the total solid region of multiple branches will be assembled by level-set functions of different Y-shape branches using R-functions.38

The coordination of the Gauss points will be substituted into the above assembled level-set equation. When the positive and negative level-set function values simultaneously exist at the Gauss points in a single element, it is considered that the element is cut by the contour and should be further subdivided. The procedure will be iteratively implemented three times for the refined “sub-elements”. The elements have been cut slightly but the level-set function values at the Gauss points are all positive or negative will be ignored in the subdivision procedure.

To implement the mechanical performance analysis, we consider the strong form of the boundary value problem for the region Ω

σ is the Cauchy stress tensor. u and u are the displacement vector and 0 the displacement vector on the boundary condition Γ . G is the body D force. n and f are the unit outward normal vector on the loading boundary Γ and the boundary traction force applied on Γ. The weak f f form for Ω can be expressed as

where ε and E denote the strain tensor and the elasticity tensor,respectively. w is a set of virtual displacements. In FCM, the unknown displacement field is expressed with high-order shape function matrix N

ω(x) is a weighting factor function, which satisfies ω(x)=0 onΓ . Based D on the above definitions, the finite element equation can be expressed as..

where B is the geometry matrix, D is the elastic matrix. The global stiffness matrix K and force vector F will be calculated by Gaussian quadrature over all elements and sub-elements. In this way, the design sensitivities of the mechanical responses respect to the feature shape and layout parameters can be immediately transferred to the design sensitivities respect to the element equivalent stiffness on the boundary, which are implemented with a boundary integral. Similar derivations can be found in some recent references.32,36

Meanwhile, the panel of the surface is modeled as shell elements that will be bonded together with the Y-shape branches using face-toface Multipoint Constraints (MPC). In this way, the direct nodeconnections are not needed, greatly simplifying the modelling procedure. Meanwhile, it should be noted that as the Y-shape branches are modelled with solid elements in the fixed background mesh, the corresponding nodes have only translational DOFs. The rotational DOFs on the skin shell elements will not be passed to the solid elements.

To verify the design procedure and the effect of the bio-inspired featuredriven topology optimization, a typical trapezoidal titanium rudder with a single shaft connecting to the body of the flight vehicle is taken as an example. The geometry and the original structural design are presented as an open rudder are shown in Fig. 7. The structural branches are distributed regularly. The basic thickness of the skin is 3mm and the core is considered as the design domain.

**Fig 7.** A trapezoidal rudder and its original design for further optimization.

**Fig 8.** Von-Mises stress distribution and global deformation of the rudder.

**Fig 9.** Optimized design obtained by topology optimization and the rebuilt structure using Y-shape branches.

The structural mass is 14.469 kg. A uniformly distributed aerodynamic pressure 0.2 Mpa is applied on the surface with a deformation as show in Fig. 8. The maximum deformation is 4.1 mm on the rear tip and the maximum Von-Mises stress is 395 Mpa near the shaft. The maximum stress on the skin is 172 Mpa. Sectional views of the stresses are provided on thickness direction of different percentages.According to the deformation pattern, the displacements on the leading edge are much smaller than those on the trailing edge, which indicates that the core structure of the existing design is much too stiff in the leading part against the trailing part. Structural optimization is therefore essential to balance the stiffness of the rudder.

To have a high-performance initial design for further bio-inspired optimization, a classical compliance-based topology optimization using standard pseudo-density variables and SIMP interpolation11 is implemented with the same loads and boundary conditions. The core of the rudder is selected as the design domain with a 35 % volume fraction constraint. The optimized design is shown in Fig. 9, where the topologically optimized design cannot be directly used. The structural topology indicates the distribution of the Y-shape branches reconstructing the core structures. The width of each branch and the thickness of the panel can be decided accordingly. This is a standard procedure to apply topology optimization into practical engineering.The rebuilt optimized design is also shown in Figure 9 that have almost identical material usage as the existing design shown in Fig. 7.

The rebuilt optimized design is further analyzed. The stress distribution and deformation is shown in Fig. 10. The maximum deformation is now 3.5 mm on the rear tip, which is 14.6 % lower than before. The maximum Von-Mises stress is 375 Mpa near the shaft, which is 5 % lower than before. The maximum stress on the skin is 151 Mpa. Sectional views of the stresses are also provided on thickness direction of different percentages. It is seen that although the reconstruction is a very subjective procedure, the topology optimization has already obviously improved the mechanical performance of the rudder.

To initiate the feature-driven optimization, the skeleton of the optimized design is drawn in Fig. 11. The red lines indicate the existing structural branches and the green lines are the branches with minimum widths which are artificially added to the structure to increasing the degrees of freedom of the optimization design spaces to obtaining better results. During the bio-inspired topology optimization, the sizes and layouts of the Y-shape branches evolve to optimal design. The iteration history is also shown in Fig. 11. Especially, some of the branches are inefficient and generally shrink and disappear in the design domain.

After 72 iterations, the optimization reached the convergence and

**Fig 10.** Von-Mises stress distribution and global deformation of the optimized rudder obtained by topology optimization.

**Fig 11.** Initial skeleton and its iteration history of the bio-inspired topology optimization.

**Fig 12.** Optimized design obtained by the bio-inspired topology optimization and the rebuilt structure using Y-shape branches.

**Fig 13.** Von-Mises stress distribution and global deformation of the optimized rudder obtained by the bio-inspired topology optimization.

the final design is shown in Fig. 12. Compared with the initial optimized design, we find that near the fixation of the rudder more material is used and a large number of branching structures are concentrated nearby. Meanwhile, the final design uses larger and longer Y-shape branches instead of some complicated multi-level Y-shape branches, which also appears in the operculum structure. We rebuilt the final structure as clear 3D configuration and analyzed it to show the effect of the optimization. The final design has almost identical material volume with the original design shown in Fig. 7 and the initial optimized design shown in Fig. 9. The maximum deformation is now only 2.9 mm on the rear tip, which is 17.1 % lower than the initial optimized design and 29.2 % lower than the original design. The maximum Von-Mises stress is 295 Mpa near the shaft, which is 21.3 % lower than the initial optimized design and 25.3 % lower than the original design. The maximum stress on the skin is now 110 Mpa. Sectional views of the stresses are also provided on thickness direction of different percentages. To have clearer comparison, these results are listed in Table 1.

Comparisons of these results are of great importance for the designers using topology optimization. Take the design that is directly rebuilt from the topologically optimized result for example, the maximum displacement and maximum stress level decreased considerably, indicating that the effect of the topology optimization is acceptable. Designers always stop here or carry out some sizing parameter optimization at best, which can only improve the structural performance in a limit scale. However, further bio-inspired topology optimization found an even better design with larger improvements in mechanical performances. By redistributing the Y-shape branches in a feature-driven way, the optimized design implies that the initial optimized design is still far from the final optimum. The reason is somehow due to the rebuilding procedure of the topologically optimized design. Designers with different experience levels always consider many manufacturing and functional requirements, consciously or unconsciously, which will significantly weaken the structural performance. Therefore, considering these rebuilt structural branches as some designable features and performing further optimization will be essential. Compared with the traditional topology optimization, the feature-driven bio-inspired topology optimization directly optimizes the detailed distribution for structural features and branches, which keeps the basic structural configuration and manufacturability and can be directly used for practical design. The optimized structural performance will thus be maintained in the final design of the engineering parts.

To verify the effect of the optimization, the three designs shown in Table 1 were fabricated with stereolithography based resin additive manufacturing, as shown in Fig. 14. The resin material is SPS6000B, which can produce strong, tough and water-resistant parts. The liquid resin is kept in constant 32 oC for its proper liquidity. The ultraviolet light with a wavelength of 355 nm is maintained 210 mw for solidification in the SPS350B stereolithography machine. The mechanical properties listed in the following table are taken from the previous work.39

Our previously implemented experiments have shown that the material is anisotropic to someextent. This is mostly due to the print direction and the procedure of solidification.40,41 As a result, all three designs were fabricated in the same direction and post-cured with ultraviolet light for 90 minutes to maintain similar material mechanical properties.

The static structure tests were performed with a 5 kN TestResourcesTM Testing machine under quasi-static conditions (see Fig. 15). An infrared digital camera (NDI OptoTRAK) was used to record the deformation of the specimen structures. The infrared displacement sensors were then attached to the fixtures of the testing machine and the four corners of the rudder. In order to simulate the aerodynamic pressure on the rudder, a water capsule was used in the experiment. The loading rate was set as 2 mm/min. All the specimen structures were tested under the same condition. The experimental force-displacement characteristics of all three rudder structures are also shown in Fig. 15, where the force is recorded as the pressure on the rudder and the displacement is the maximum one on the rudder corners.

Among the three designs, it can be seen that the forcedisplacement curve of the structure obtained from the bio-inspired feature-driven topology optimization has much larger slope than the other two curves. It indicates that the structure has the best stiffness performance, followed by the structure directly reconstructed from the classic topology optimization and then the original design. Hence, the experimental results are in good agreement with the results of numerical simulation shown in Table 1.

**Fig 14.** Three different designs fabricated with stereolithography based resin additive manufacturing.

**Fig 15.** Setup of the static loading experiment and the load-displacement curve of the rudders.

In this work, we focused on flight vehicle rudder structure design inspired by a fish operculum bone structure. Its typical configuration that is composed of Y-shape branches of different sizes and skeleton contours gives it excellent performance transferring surface loads to the joint position. It is also found that the distribution of the branches on the operculum bone well corresponds to the topologically optimized results. Therefore, we developed a bio-inspired topology optimization approach for flight vehicle rudder and other thin-walled structures design. The Yshape branches are considered as basic features and parameters such as their distributions on the panel, skeleton shapes and their sizes in different positions are introduced as design variables simultaneously. A feature-driven optimization scheme is thus established minimizing the structural strain energy and constraining maximum deformation and Von-Mises stresses to guarantee a stiffness and strength design. Typically, to improve the quality of the total optimization procedure, the bio-inspired feature-driven optimization scheme will start from an initial design that is obtained and reconstructed from a classic topology optimization.

A practical example was presented to demonstrate the validity and advantages of the proposed method. Compared with the original design and the one obtained directly from classic topology optimization, the performance of the bio-inspired design was significantly improved, which implies that the optimization shall not stop with the results reconstructed directly from topology optimization. Further featuredriven optimization will continue improving the structural performance on the one hand, and maintain the basic structural configuration and manufacturability on the other hand. In addition, three different designs of the rudder were fabricated by stereolithography based resin additive manufacturing and then the mechanical tests were performed under quasi-static conditions. The experimental results also confirm that the optimized design obtained by bio-inspired feature-driven optimization can improve the stiffness and strength significantly compared with the direct topology optimization design.

**Acknowledgements**

This work is supported by National Key Research and Development Program (2017YFB1102800), NSFC for Excellent Young Scholars (11722219) and Key Project of NSFC (51790171, 5171101743, 51735005).

1. C. Y. Niu, Airframe structural design: practical design information and data on flight vehicle structures, Conmilit Press Limited, 1999.

2. R. Vitali, O. Park, R. T. Haftka, B. V. Sankar, C. A. Rose, J. Aircraft., 2002, 39, 158-166.

3. S. Chintapalli, M. S. A. Elsayed, R. Sedaghati, M. Abdo, Aerosp. Sci. Technol., 2010, 14, 188-198.

4. L. U. Hansen, P. Horst, Comput. Sturct., 2008, 86, 104-118.

5. S. Grihon, D. Bassir, L. Krog, Ijsmdo., 2009, 3, 432-442.

6. B. Wang, P. Hao, G. Li, K. Tian, K. Du, X. Wang, X. Zhang, X. Tang, Struct. Multidiscip. O., 2014, 50, 313-327.

7. J. K. Allen, S. Azarm, T. W. Simpson, J. Mech. Design., 2011, 133, 100301.

8. J. H. Zhu, W. H. Zhang, L. Xia, Arch. Comput. Method. E., 2016, 23, 595-622.

9. K. T. Cheng, N. Olhoff, Int. J. Solids. Struct., 1981, 17, 305-323.

10. M. P. B. Xf, E. Kikuchi, Comput. Method. Appl. M., 1988, 71, 197- 224.

11. O. Sigmund, K. Maute, Struct. Multidiscip. O., 2013, 48, 1031- 1055.

12. J. D. Deaton, R. V. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post 2000, Springer-Verlag New York, Inc., 2014.

13. D. Walker, D. Liu, 56th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2015.

14. Y. B. Zhao, W. J. Guo, S. H. Duan, L. G. Xing, Ijsmdo., 2017, 8, A5.

15. N. Aage, E. Andreassen, B. S. Lazarov, O. Sigmund, Nature, 2017, 550, 84-86.

16. J. Hou, J. Zhu, F. He, W. Zhang, W. Guo, Chinese. J. Aeronaut., 2017, 30, 1441-1450.

17. B. K. Stanford, P. D. Dunning, J Aircraft, 2015, 52, 1298-1311.

18. J. H. Zhu, J. Hou, W. H. Zhang, Y. Li, Struct. Multidiscip. O., 2014, 50, 561-571.

19. J. Luo, H. C. Gea, Struct. Optim., 1998, 16, 280-288.

20. H. C. Gea, J. Luo, Mech. Struct. Mach., 1999, 27, 275-292.

21. D. Wang, M. M. Abdalla, W. Zhang, Compos. Struct., 2017, 159, 656-666.

22. D. Wang, M. M. Abdalla, W. Zhang, Compos. Struct., 2018, 193, 224-236.

23. G. C. Young, G. R. Zhang, Palaeontology, 1992, 35, 443-464.

24. C. B. Kimmel, W. E. Aguirre, B. Ullmann, M. Currey, W. A. Cresko, Behaviour, 2008, 145, 669-691.

25. J. X. Zhang, W. Bo, N. Fei, G. D. Cheng, Chin. J. Comput. Mech., 2014, 31, 141-148.

26. Y. Cao, J. Zhu, W. Guo, Y. Li, W. Zhang, Sci. Sinica, 2018, 48, 014606.

27. T. Gao, W. H. Zhang, J. H. Zhu, Y. J. Xu, D. H. Bassir, Finite. Elem. Anal. Des., 2008, 44, 805-813.

28. A. Iga, S. Nishiwaki, K. Izui, M. Yoshimura, Int. J. Heat. Mass. Tran., 2009, 52, 2721-2732.

29. J. R. R. A. Martins, N. M. K. Poon, VI World Congress on Structural and Multidisciplinary Optimization WCSMO6, Rio de Janeiro, Brasil. Citeseer, 2005.

30. N. M. K. Poon, J. R. R. A. Martins, Struct. Multidiscip. O., 2007, 34, 61-73.

31. H. H. Gao, J. H. Zhu, W. H. Zhang, Y. Zhou, Comput. Method. Appl. M., 2015, 289, 387-408.

32. Y. Zhou, W. H. Zhang, J. H. Zhu, Z. Xu, Comput. Method. Appl. M., 2016, 310, 1-32.

33. X. Guo, W. Zhang, W. Zhong, Comput. Method. Appl. M., 2014, 272, 354-378.

34. P. Wei, M. Y. Wang, X. Xing, Comput. Aided. Design., 2010, 42, 708-719.

35. A. Düster, J. Parvizian, Z. Yang, E. Rank, Comput. Method. Appl. M., 2008, 197, 3768-3782.

36. S. Y. Cai, W. H. Zhang, J. H. Zhu, T. Gao, Comput. Method. Appl. M., 2014, 278, 361-387.

37. L. Xia, J. H. Zhu, W. H. Zhang, Comput. Method. Appl. M., 2012, 241-244, 142-154.

38. V. Shapiro, Technical Report, Cornell University, 1991.

39. C. Wang, J. H. Zhu, W. H. Zhang, S. Y. Li, J. Kong, Struct. Multidiscip. O., 2018, 58, 35-50.

40. K. K. Yang, J. H. Zhu, C. Wang, D. S. Jia, L. L. Song, W. H. Zhang, Comput. Mech., 2018, 61, 581-598.

41. T. Liu, S. Guessasma, J. Zhu, W. Zhang, H. Nouri, S. Belhabib, J. Mater. Process. Tech., 2017, 251, 37-46.