Bionic Stiffener Layout Optimization with a Flexible Plate in Solar-Powered UAV Surface Structure Design

A cellular-based evolutionary topology optimization scheme over a small curvature big contour wing surface is proposed for the design of an ultralight surface structure. Using this method, a ground-structure technique is first applied to obtain homogeneous mesh generation with a predefined weight value over the design domain. Secondly, the stiffener path’s description is guided by a modified map L system topology method that simulates the growth of the bionic branch, and the structural components are obtained by the specified searching method according to weights of the previous mesh vertexes. Thirdly, an optimal curved stiffener layout is achieved using an agent-based algorithm to create individual instances of designs based on a small number of input parameters. These parameters can then be controlled by a genetic algorithm to optimize the final design according to goals like minimizing weight and structural weakness. A comparison is implemented for long-span panel stiffener layout generation between an initial straight case and a bionic optimal case via our method, thereby indicating the significant improvement of the buckling loads by steering the stiffener’s path. Finally, this bionic method is applied to the wing box structure design and achieves remarkable weight loss at last.


Introduction
Solar-electrically powered fixed-wing unmanned aerial vehicles (UAVs), regarded as a substitute for satellites, have wide applications in communication, disaster relief, supervisory control, and so on. In order to obtain a larger lift-to-drag ratio and solar panel area, this vehicle usually has a high-aspect ratio wing, which requires a large size and very flexible wing-structure. Because of their strict weight limitations, these UAVs often have extreme structural optimization requirements to balance their deformation and material distribution. With the development of solar cell packaging technology, shown in Figure 1, these solar cell modules have become much lighter and more flexible than before, which also means that these structures may fail before the wing achieves maximum tip displacement owing to local buckling on its compressive side. Thus, it is challenging to both minimize the structure's weight and improve the performance of the structure.
Laminated composite structures with unidirectional glass or carbon fiber reinforced polymer (GFRP/CFRP) materials are usually used as a substrate and are connected together by a resin with solar cell modules to improve their strength and stiffness characteristics [1][2][3]. Generally, this laminated structure with a large thickness may cause the weight to exceed its limits, while a thinner thickness may reduce its stiffness. Consequently, a large number of cutouts exist in this thin-wall structure as a direct method to save the weight. However, the presence of cutouts will destroy the integrity of the structure, reducing buckling resistance and also possibly generating unpredictable flaws. In order to obtain a lighter weight structure, the traditional optimization method is to optimize the laminate stacking sequence. As the development of automated fiber placement techniques makes curved fiber paths affordable, the concept of variable stiffness laminates with curved fiber paths is proposed. With the development of manufacturing technologies, more efficient structural forms have become available. These new structures offer tremendous opportunities for combination with optimization technologies to obtain satisfactory structural weight reduction.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 17 thickness may reduce its stiffness. Consequently, a large number of cutouts exist in this thin-wall structure as a direct method to save the weight. However, the presence of cutouts will destroy the integrity of the structure, reducing buckling resistance and also possibly generating unpredictable flaws. In order to obtain a lighter weight structure, the traditional optimization method is to optimize the laminate stacking sequence. As the development of automated fiber placement techniques makes curved fiber paths affordable, the concept of variable stiffness laminates with curved fiber paths is proposed. With the development of manufacturing technologies, more efficient structural forms have become available. These new structures offer tremendous opportunities for combination with optimization technologies to obtain satisfactory structural weight reduction. As illustrated in Figure 2, grid-stiffened structures or shell structures preserve a high strength because of their exceptional mechanical properties along the fiber direction and are now widely used in both aerospace engineering and ground transportation vehicles. On the basis of the concept of variable stiffness laminates with curved fiber paths proposed by Gürdal [4], a new idea for arbitrarily distributed curved stiffener layouts has inspired a new wave of research. Dan Wang made full use of these streamlined characteristic to design a fiber path, and these non-uniform curved grid-stiffened composite (NCGC) structures have achieved a remarkable increase in their buckling load [5,6]. Kapania engaged in extensive research on the shape optimization of curved stiffeners for composite structures [7], while Mohammad Rouhi proposed a variable stiffness composite cylinder by fiber steering to maximum the buckling load under pure bending [8]. Beyond the shape or size optimization of the stiffeners or varying the path of the stiffeners' distribution, there are several novel features in this paper's surface design. Firstly, this design offers a more flexible surface. Secondly, the grid-like structure is developed, with more grids added in the high stress area to improve safety margins. Thirdly, this design features a non-uniform load distribution. These particularities have yielded a new concept distinct from current optimization methods and also offer an opportunity to determine new structural types with high efficiency. As illustrated in Figure 2, grid-stiffened structures or shell structures preserve a high strength because of their exceptional mechanical properties along the fiber direction and are now widely used in both aerospace engineering and ground transportation vehicles. On the basis of the concept of variable stiffness laminates with curved fiber paths proposed by Gürdal [4], a new idea for arbitrarily distributed curved stiffener layouts has inspired a new wave of research. Dan Wang made full use of these streamlined characteristic to design a fiber path, and these non-uniform curved grid-stiffened composite (NCGC) structures have achieved a remarkable increase in their buckling load [5,6]. Kapania engaged in extensive research on the shape optimization of curved stiffeners for composite structures [7], while Mohammad Rouhi proposed a variable stiffness composite cylinder by fiber steering to maximum the buckling load under pure bending [8].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 17 thickness may reduce its stiffness. Consequently, a large number of cutouts exist in this thin-wall structure as a direct method to save the weight. However, the presence of cutouts will destroy the integrity of the structure, reducing buckling resistance and also possibly generating unpredictable flaws. In order to obtain a lighter weight structure, the traditional optimization method is to optimize the laminate stacking sequence. As the development of automated fiber placement techniques makes curved fiber paths affordable, the concept of variable stiffness laminates with curved fiber paths is proposed. With the development of manufacturing technologies, more efficient structural forms have become available. These new structures offer tremendous opportunities for combination with optimization technologies to obtain satisfactory structural weight reduction. As illustrated in Figure 2, grid-stiffened structures or shell structures preserve a high strength because of their exceptional mechanical properties along the fiber direction and are now widely used in both aerospace engineering and ground transportation vehicles. On the basis of the concept of variable stiffness laminates with curved fiber paths proposed by Gürdal [4], a new idea for arbitrarily distributed curved stiffener layouts has inspired a new wave of research. Dan Wang made full use of these streamlined characteristic to design a fiber path, and these non-uniform curved grid-stiffened composite (NCGC) structures have achieved a remarkable increase in their buckling load [5,6]. Kapania engaged in extensive research on the shape optimization of curved stiffeners for composite structures [7], while Mohammad Rouhi proposed a variable stiffness composite cylinder by fiber steering to maximum the buckling load under pure bending [8]. Beyond the shape or size optimization of the stiffeners or varying the path of the stiffeners' distribution, there are several novel features in this paper's surface design. Firstly, this design offers a more flexible surface. Secondly, the grid-like structure is developed, with more grids added in the high stress area to improve safety margins. Thirdly, this design features a non-uniform load distribution. These particularities have yielded a new concept distinct from current optimization methods and also offer an opportunity to determine new structural types with high efficiency. Beyond the shape or size optimization of the stiffeners or varying the path of the stiffeners' distribution, there are several novel features in this paper's surface design. Firstly, this design offers a more flexible surface. Secondly, the grid-like structure is developed, with more grids added in the high stress area to improve safety margins. Thirdly, this design features a non-uniform load distribution. These particularities have yielded a new concept distinct from current optimization methods and also offer an opportunity to determine new structural types with high efficiency. The remainder of the paper is organized as follows. In Section 2, a bionic concept for layout generation is proposed. Next, in Section 3, the details of the bionic stiffener's growth process are presented. The optimization settings and discussions of efficiency, followed by FE (Finite Element) analysis, are listed in Section 4. In Section 5, a wing structure established by this article's method is discussed. In Section 6, we conclude this paper with some remarks.

Bionic Concept of the Grid Generation Method
What seems challenging to us is very common in some natural structures, such as insect wings, plant leaves, and diatom shells, as shown in Figure 3. These thin-wall structures simultaneously use multiple intersected curved stiffeners to form their surfaces and achieve great enough strength to protect themselves from outside damage.
What seems challenging to us is very common in some natural structures, such as insect wings, plant leaves, and diatom shells, as shown in Figure 3. These thin-wall structures simultaneously use multiple intersected curved stiffeners to form their surfaces and achieve great enough strength to protect themselves from outside damage.
Many designers are inspired by these natural structures and have applied bionic design methods to obtain optimal solutions. Jihong Zhu demonstrated a shape preserving topology optimization design approach for suppressing the warping deformation of local structural domains [9]. Dimcic Milos applied a relaxed method based on the force density method to a non-uniform Voronoi structure and obtained a foam-like structure through the genetic algorithm [10]. Danil Nagyand and this APworks group combined a bottom-up growth strategy based on slime mould's behavior in nature using a top-down genetic algorithm strategy and by minimizing weight and structural weakness to obtain the final design [11,12]. Christian Hammand and his group have engaged in extensive research on plankton shells and presented an ELiSE (evolutionary light structure engineering) concept that uses pre-optimized lightweight structures to widen the design space through the development of various new lightweight solutions to determine the best type of structure [13,14]. These studies have established practical indirect strategies to simulate the characteristics of natural structures by changing their parameters to produce complex design solutions that are not only high performing, but also novel and unexpected. To take advantage of nature's structure and simulate growth behavior, a two-step stiffener growth method is presented to combine the biomimetic domain partition method with the structural topology method and insert an evolution algorithm into the design process to obtain an optimal result. As shown in Figure 4, a rule-based algorithm is used to obtain various stiffener layouts for a given wing surface geometry. First, a predicted topology method based on the ground structure (GS) is applied to obtain the discrete structured data. Then, a modified Lindenmeyer systems (L-systems) algorithm based on the homogenized structural data is executed to create the actual structure. A detailed structural model is established by parametric modeling, and an FEM simulation is implemented to obtain the required judgment data for the optimization algorithm. Many designers are inspired by these natural structures and have applied bionic design methods to obtain optimal solutions. Jihong Zhu demonstrated a shape preserving topology optimization design approach for suppressing the warping deformation of local structural domains [9]. Dimcic Milos applied a relaxed method based on the force density method to a non-uniform Voronoi structure and obtained a foam-like structure through the genetic algorithm [10]. Danil Nagyand and this APworks group combined a bottom-up growth strategy based on slime mould's behavior in nature using a top-down genetic algorithm strategy and by minimizing weight and structural weakness to obtain the final design [11,12]. Christian Hammand and his group have engaged in extensive research on plankton shells and presented an ELiSE (evolutionary light structure engineering) concept that uses pre-optimized lightweight structures to widen the design space through the development of various new lightweight solutions to determine the best type of structure [13,14]. These studies have established practical indirect strategies to simulate the characteristics of natural structures by changing their parameters to produce complex design solutions that are not only high performing, but also novel and unexpected.
To take advantage of nature's structure and simulate growth behavior, a two-step stiffener growth method is presented to combine the biomimetic domain partition method with the structural topology method and insert an evolution algorithm into the design process to obtain an optimal result. As shown in Figure 4, a rule-based algorithm is used to obtain various stiffener layouts for a given wing surface geometry. First, a predicted topology method based on the ground structure (GS) is applied to obtain the discrete structured data. Then, a modified Lindenmeyer systems (L-systems) algorithm based on the homogenized structural data is executed to create the actual structure. A detailed structural model is established by parametric modeling, and an FEM simulation is implemented to obtain the required judgment data for the optimization algorithm.

Modeling Framework
Here, we provide a brief introduction to the ground structure method and the L-systems. A detailed process introduction is also presented for the two-step stiffeners growth method. There are three basic points in this strategy to form a bionic structure layout; all points are separately arranged into two-steps: (1) discretization, (2) homogenization, and (3) partitioning. With these factors, the parameters are extracted and used to control the rules of the intermediate algorithm, which is then executed to create the actual structure.

Ground Structure Method
The ground-structure technique, which was first introduced in the 1960s by Dorn, is usually used for the truss topology design problem. In this approach, a large number of potential nodes and an even larger number of potential bars are distributed over a design domain. With a continuum of geometrically restricted conditions, this sizing reformulation ultimately generates the truss-layout by removing unnecessary members [15]. A stable truss system that is sufficient to prevent rigid body motion should satisfy the force equilibrium equations to obtain the minimum volume. The formulation is where V is the the truss' volume; and , , and are the cross-sectional area, length, and stress, respectively, of the ℎ truss member (or all members in the truss). is internal bar forces. The parameters and are the numbers of nodes and components with supports, respectively, and are the free nodal components with = n − for n = 2 in a two-dimensional ground structure.
is the nodal equilibrium matrix of size × , is the force's node, and its

Modeling Framework
Here, we provide a brief introduction to the ground structure method and the L-systems. A detailed process introduction is also presented for the two-step stiffeners growth method. There are three basic points in this strategy to form a bionic structure layout; all points are separately arranged into two-steps: (1) discretization, (2) homogenization, and (3) partitioning. With these factors, the parameters are extracted and used to control the rules of the intermediate algorithm, which is then executed to create the actual structure.

Ground Structure Method
The ground-structure technique, which was first introduced in the 1960s by Dorn, is usually used for the truss topology design problem. In this approach, a large number of potential nodes and an even larger number of potential bars are distributed over a design domain. With a continuum of geometrically restricted conditions, this sizing reformulation ultimately generates the truss-layout by removing unnecessary members [15]. A stable truss system that is sufficient to prevent rigid body motion should satisfy the force equilibrium equations to obtain the minimum volume. The formulation is where V is the the truss' volume; and a i , l i , and σ i are the cross-sectional area, length, and stress, respectively, of the ith truss member (or all N b members in the truss). n is internal bar forces. The parameters N n and N sup are the numbers of nodes and components with supports, respectively, and N do f are the free nodal components with N do f = n N n − N sup for n = 2 in a two-dimensional ground structure. B T is the nodal equilibrium matrix of size N do f × N b , f is the force's node, and its number is N do f ; n is a vector with an internal (axial) force for all members in the ground structure. Stress limits in tension are σ T > 0 and compression is σ C < 0.
To convert the constraint into its equalities, the slack variables s + and s − are used, and the optimization problem in (1) becomes a linear programming problem, as follows: The member experiences tension if s + i > 0 and compression if s + i > 0. Using matrix notation, (2) can be rewritten as Equation (3) and can be solved more efficiently using the interior-point algorithm [16,17]: There are three major reasons to choose the ground-structure technique as the first topology method to reflect the layout pattern of the given surface domain: (i) the use of a plate-like design domain for the wing surface; (ii) the length of this columnar structure is much greater than its dimensions (diameter or side lengths), and its prime characteristics are along the carbon fiber, which acts like a truss structure; (iii) evenly distributed grid vertexes, which can discretize the whole structure's information and function as knots along the detail stiffeners' path.

Discrete Strategy
In this strategy, an amount of discrete points among the design domain is needed to contain the structural data after the topology evaluation step under the GS method, as the GS is a polygon mesh-based method, which means a ruled partition is used for the given design domain. The usual used grid mesh for the GS is rectangular and triangular, with polygon vertexes that are constant in the domain. Centroidal Voronoi tessellation (CVT) is a special condition of the Voronoi diagram, in which the generating point pi coincides with the centroid pic of the corresponding region (shown in Figure 5) [18]. This has wide applications in data compression, cell biology, the territorial behavior of animals, and the optimal allocation of resources. This special convex polygon has a random total edge range from 4 to 6, and their corresponding vertexes enlarge the search range of the truss orientation in the GS, which is advantageous for the searching strategy during the stiffener growth.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 5 of 17 number is ; is a vector with an internal (axial) force for all members in the ground structure. Stress limits in tension are > 0 and compression is < 0. To convert the constraint into its equalities, the slack variables + and − are used, and the optimization problem in (1) where = + / + − / and = + − − , = / . Only one of + and − is non-zero. The member experiences tension if + > 0 and compression if + > 0. Using matrix notation, (2) can be rewritten as Equation (3) and can be solved more efficiently using the interior-point algorithm [16][17]: There are three major reasons to choose the ground-structure technique as the first topology method to reflect the layout pattern of the given surface domain: (i) the use of a plate-like design domain for the wing surface; (ii) the length of this columnar structure is much greater than its dimensions (diameter or side lengths), and its prime characteristics are along the carbon fiber, which acts like a truss structure; (iii) evenly distributed grid vertexes, which can discretize the whole structure's information and function as knots along the detail stiffeners' path.

Discrete Strategy
In this strategy, an amount of discrete points among the design domain is needed to contain the structural data after the topology evaluation step under the GS method, as the GS is a polygon meshbased method, which means a ruled partition is used for the given design domain. The usual used grid mesh for the GS is rectangular and triangular, with polygon vertexes that are constant in the domain. Centroidal Voronoi tessellation (CVT) is a special condition of the Voronoi diagram, in which the generating point pi coincides with the centroid pic of the corresponding region (shown in Figure 5) [18]. This has wide applications in data compression, cell biology, the territorial behavior of animals, and the optimal allocation of resources. This special convex polygon has a random total edge range from 4 to 6, and their corresponding vertexes enlarge the search range of the truss orientation in the GS, which is advantageous for the searching strategy during the stiffener growth.  In this strategy, in order to obtain a 3D geometric mesh, a non uniform rational basis spline (NURBS) method is used by Rhino (3D model built software) to transform the parameters from a two-dimension plane (x, y) to the referenced three-dimensional surface (u, v). The basic mesh for the wing surface of the GRAND method is shown in Figure 6. After establishing the plane grid's relationship and setting the boundary and load conditions for the mesh, we use the optimization function defined in Equation (3) to obtain the bar connections and the bar size, respectively.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 17 In this strategy, in order to obtain a 3D geometric mesh, a non uniform rational basis spline (NURBS) method is used by Rhino (3D model built software) to transform the parameters from a two-dimension plane (x, y) to the referenced three-dimensional surface (u, v). The basic mesh for the wing surface of the GRAND method is shown in Figure 6. After establishing the plane grid's relationship and setting the boundary and load conditions for the mesh, we use the optimization function defined in Equation (3) to obtain the bar connections and the bar size, respectively.

Homogenization Strategy
In this strategy, topological information obtained from the discrete strategy will be homogenized. As depicted in Figure 7, The GS computed a layout configuration that also reflects the material's distribution over the design domain. For example, in the pre-defined points of the design domain, a greater number of bars at the same point means more materials in this area, which also indicates a high stress level. Therefore, this homogenization method is a size parameter homogenized process for arranging the bar members' materials to the corresponding point. As one location may belong to two or more bars' ends, and the locations can also be influenced by the surrounding bars, three steps are executed to homogenize the truss connections to the mass data. Figure 8 shows a separate process for homogenizing the bar member's materials to the corresponding points. The corresponding points are selected by the following steps: (1) The searching domain is established among the two ends of the bar.
(2) Calculating the distance from the corresponding point to the bar member.
(3) Allocating the size parameter associated with the respective distance. (4) Data summation and normalization processing for every point, : where is the maximum computed weight data, and W is the weight from the ith corresponding bar. After the above four steps, the topology data will change from the bar size data to the stress-weight data of the points (seeds) in the domain.

Homogenization Strategy
In this strategy, topological information obtained from the discrete strategy will be homogenized. As depicted in Figure 7, The GS computed a layout configuration that also reflects the material's distribution over the design domain. For example, in the pre-defined points of the design domain, a greater number of bars at the same point means more materials in this area, which also indicates a high stress level. Therefore, this homogenization method is a size parameter homogenized process for arranging the bar members' materials to the corresponding point. In this strategy, in order to obtain a 3D geometric mesh, a non uniform rational basis spline (NURBS) method is used by Rhino (3D model built software) to transform the parameters from a two-dimension plane (x, y) to the referenced three-dimensional surface (u, v). The basic mesh for the wing surface of the GRAND method is shown in Figure 6. After establishing the plane grid's relationship and setting the boundary and load conditions for the mesh, we use the optimization function defined in Equation (3) to obtain the bar connections and the bar size, respectively.

Homogenization Strategy
In this strategy, topological information obtained from the discrete strategy will be homogenized. As depicted in Figure 7, The GS computed a layout configuration that also reflects the material's distribution over the design domain. For example, in the pre-defined points of the design domain, a greater number of bars at the same point means more materials in this area, which also indicates a high stress level. Therefore, this homogenization method is a size parameter homogenized process for arranging the bar members' materials to the corresponding point. As one location may belong to two or more bars' ends, and the locations can also be influenced by the surrounding bars, three steps are executed to homogenize the truss connections to the mass data. Figure 8 shows a separate process for homogenizing the bar member's materials to the corresponding points. The corresponding points are selected by the following steps: (1) The searching domain is established among the two ends of the bar.
(2) Calculating the distance from the corresponding point to the bar member.
(3) Allocating the size parameter associated with the respective distance. (4) Data summation and normalization processing for every point, : where is the maximum computed weight data, and W is the weight from the ith corresponding bar. After the above four steps, the topology data will change from the bar size data to the stress-weight data of the points (seeds) in the domain. As one location may belong to two or more bars' ends, and the locations can also be influenced by the surrounding bars, three steps are executed to homogenize the truss connections to the mass data. Figure 8 shows a separate process for homogenizing the bar member's materials to the corresponding points. The corresponding points are selected by the following steps: (1) The searching domain is established among the two ends of the bar.
(2) Calculating the distance from the corresponding point to the bar member.
(3) Allocating the size parameter associated with the respective distance. (4) Data summation and normalization processing for every point, w e : where w e max is the maximum computed weight data, and W i is the weight from the ith corresponding bar. After the above four steps, the topology data will change from the bar size data to the stress-weight data of the points (seeds) in the domain. Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 17 Figure 8. An allocation method to homogenize the size data from the bar member to the selected node.

Map L Systems
To generate an entire grid structure with no disconnection, an intermediate system is needed to guide the partition under a certain rule. Thus, a binary propagating map Lindenmeyer system, whose grid follows certain rules, is used here [19]. In this system, the domain partitions and cell edges are generated by a set of rules. One example of the production rules is presented in Equation (5), and the derivation process of the first four steps is given in Figure 9: where  

,,
A B x  are letters that define the class of the edges. Detailed explanations of the rules to this equation can be found in the literature [20][21][22]. The process may terminate after each cell is less than the required cell area, after a defined maximum number of steps, or after certain production rules terminate on their own volition. The layout of the final cell is finally transformed into the structural-design domain. In our method, in order to utilize the information from the GS step, a modified cellular-division methodology is proposed and detailed in the partition strategy.

Partition Strategy
In this strategy, a modified L-system based on GS is presented for a non-uniform curved stiffener layout design. After the homogenization strategy, every location contains the weight data. To enact the L-system in our design model, we first specified two types of locations, the 'edge seeds' and 'body seeds', respectively. As shown in Figure 10, the 'edge seeds' consist of the boundary locations that would change the initial weight value, and the 'body seeds' comprise the locations in the domain with the weight value unchanged from the homogenization strategy. Before binary partition, the longer coupled edges in triangle or quadrangle are analyzed to determine the partition objects. Then, a node is chosen, and a stiffener growth process is executed to imitate the logic of the L-system. Detailed implementations for the initial binary partition domain are shown in Figure 11. It should be noted that a decay operator is needed to decay the weight of the upper level's chosen node weight.

Map L Systems
To generate an entire grid structure with no disconnection, an intermediate system is needed to guide the partition under a certain rule. Thus, a binary propagating map Lindenmeyer system, whose grid follows certain rules, is used here [19]. In this system, the domain partitions and cell edges are generated by a set of rules. One example of the production rules is presented in Equation (5), and the derivation process of the first four steps is given in Figure 9: where = {A, B, x} are letters that define the class of the edges. Detailed explanations of the rules to this equation can be found in the literature [20][21][22]. The process may terminate after each cell is less than the required cell area, after a defined maximum number of steps, or after certain production rules terminate on their own volition. The layout of the final cell is finally transformed into the structural-design domain. In our method, in order to utilize the information from the GS step, a modified cellular-division methodology is proposed and detailed in the partition strategy.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 17 Figure 8. An allocation method to homogenize the size data from the bar member to the selected node.

Map L Systems
To generate an entire grid structure with no disconnection, an intermediate system is needed to guide the partition under a certain rule. Thus, a binary propagating map Lindenmeyer system, whose grid follows certain rules, is used here [19]. In this system, the domain partitions and cell edges are generated by a set of rules. One example of the production rules is presented in Equation (5), and the derivation process of the first four steps is given in Figure 9: where  

,,
A B x  are letters that define the class of the edges. Detailed explanations of the rules to this equation can be found in the literature [20][21][22]. The process may terminate after each cell is less than the required cell area, after a defined maximum number of steps, or after certain production rules terminate on their own volition. The layout of the final cell is finally transformed into the structural-design domain. In our method, in order to utilize the information from the GS step, a modified cellular-division methodology is proposed and detailed in the partition strategy.

Partition Strategy
In this strategy, a modified L-system based on GS is presented for a non-uniform curved stiffener layout design. After the homogenization strategy, every location contains the weight data. To enact the L-system in our design model, we first specified two types of locations, the 'edge seeds' and 'body seeds', respectively. As shown in Figure 10, the 'edge seeds' consist of the boundary locations that would change the initial weight value, and the 'body seeds' comprise the locations in the domain with the weight value unchanged from the homogenization strategy. Before binary partition, the longer coupled edges in triangle or quadrangle are analyzed to determine the partition objects. Then, a node is chosen, and a stiffener growth process is executed to imitate the logic of the L-system. Detailed implementations for the initial binary partition domain are shown in Figure 11. It should be noted that a decay operator is needed to decay the weight of the upper level's chosen node weight.

Partition Strategy
In this strategy, a modified L-system based on GS is presented for a non-uniform curved stiffener layout design. After the homogenization strategy, every location contains the weight data. To enact the L-system in our design model, we first specified two types of locations, the 'edge seeds' and 'body seeds', respectively. As shown in Figure 10, the 'edge seeds' consist of the boundary locations that would change the initial weight value, and the 'body seeds' comprise the locations in the domain with the weight value unchanged from the homogenization strategy. Before binary partition, the longer coupled edges in triangle or quadrangle are analyzed to determine the partition objects. Then, a node is chosen, and a stiffener growth process is executed to imitate the logic of the L-system. Detailed implementations for the initial binary partition domain are shown in Figure 11. It should be noted that a decay operator is needed to decay the weight of the upper level's chosen node weight. In this way, the chosen node will not be the second selected in the next level of the L-system. Lastly, the binary partition domain will not be split after the cell area is less than a certain value.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 17 In this way, the chosen node will not be the second selected in the next level of the L-system. Lastly, the binary partition domain will not be split after the cell area is less than a certain value.

Model Extraction and Finite Element Analysis Preparation
In order to illustrate the potential of this bionic stiffener growing method for the buckling design of plate-like wing surface structures, we use long-span panel buckling optimization as the verification. The shell model is made of a honeycomb sandwich structure with a length of 1.5 m and a width of 0.2643 m. The up and down sandwich composite's ply thickness is 0.21 mm, and the honeycomb's core thickness is 1.58 mm; the total thickness is 2.0 mm. The two loading cases of uniaxial compression and biaxial compression are illustrated in Figure 12   In this way, the chosen node will not be the second selected in the next level of the L-system. Lastly, the binary partition domain will not be split after the cell area is less than a certain value.

Model Extraction and Finite Element Analysis Preparation
In order to illustrate the potential of this bionic stiffener growing method for the buckling design of plate-like wing surface structures, we use long-span panel buckling optimization as the verification. The shell model is made of a honeycomb sandwich structure with a length of 1.5 m and a width of 0.2643 m. The up and down sandwich composite's ply thickness is 0.21 mm, and the honeycomb's core thickness is 1.58 mm; the total thickness is 2.0 mm. The two loading cases of uniaxial compression and biaxial compression are illustrated in Figure 12

Model Extraction and Finite Element Analysis Preparation
In order to illustrate the potential of this bionic stiffener growing method for the buckling design of plate-like wing surface structures, we use long-span panel buckling optimization as the verification. The shell model is made of a honeycomb sandwich structure with a length of 1.5 m and a width of 0.2643 m. The up and down sandwich composite's ply thickness is 0.21 mm, and the honeycomb's core thickness is 1.58 mm; the total thickness is 2.0 mm. The two loading cases of uniaxial compression and biaxial compression are illustrated in Figure 12 [23].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 17 In this way, the chosen node will not be the second selected in the next level of the L-system. Lastly, the binary partition domain will not be split after the cell area is less than a certain value.

Model Extraction and Finite Element Analysis Preparation
In order to illustrate the potential of this bionic stiffener growing method for the buckling design of plate-like wing surface structures, we use long-span panel buckling optimization as the verification. The shell model is made of a honeycomb sandwich structure with a length of 1.5 m and a width of 0.2643 m. The up and down sandwich composite's ply thickness is 0.21 mm, and the honeycomb's core thickness is 1.58 mm; the total thickness is 2.0 mm. The two loading cases of uniaxial compression and biaxial compression are illustrated in Figure 12 Table 1, and 45 • and −45 • intersected straight stiffeners are used as the initial stiffener layout (illustrated in Figure 13). These two structure types (initial straight and bionic optimal) have the same skin layout and the same stiffener height. Each beam-like member of the curved stiffener has a uniform cross section. The stiffener's elastic modulus is found from the Halpin-Tsai semi-empirical relation because of its anisotropic and composite characteristics [24]. The engineering constants used to define the beam element are listed in Table 1, and 45° and −45° intersected straight stiffeners are used as the initial stiffener layout (illustrated in Figure 13). These two structure types (initial straight and bionic optimal) have the same skin layout and the same stiffener height.

Research Model Parameters
The final plate's structural design is defined by the boundary of the panel plus its partition edges. The parameters of this model are the weights of the 'boundary seeds' ( ) that are given as real values in the domain [0, 1], as well as the weights of the 'body seeds' ( ) that are calculated from the homogenization strategy, the judgement of the minimum cell area (area), the highest level number of the L-system (N), and the decay parameter (d). Here, we chose as the changed parameters; the structure will be generated under the defined rules. Because all the input parameters are continuous, the genetic algorithm (GA) is able to 'learn' how to work with the growth behavior and tune it to create better performing designs over time.

Optimization Model and Option Settings
This behavioral generative geometry model can create numerous layout patterns for a wing surface structure based on a relatively small set of input parameters (87 boundary seeds). To find high-performing layout patterns, a series of judgments are needed to the algorithm. In this paper, static finite element analysis (FEA) is used to simulate the performance of each design under the given loading conditions. This analysis gives us a set of metrics that we can use to establish the objectives and constraints of our optimization problem:  The total weight of the whole model must be minimum, with the first buckling loads confined to a certain range.  The first buckling loads, which are the objective of this design, should be maximum. These two optimization problems are presented as follows.

Research Model Parameters
The final plate's structural design is defined by the boundary of the panel plus its partition edges. The parameters of this model are the weights of the 'boundary seeds' (w bc ) that are given as real values in the domain [0, 1], as well as the weights of the 'body seeds' (w initial ) that are calculated from the homogenization strategy, the judgement of the minimum cell area (area), the highest level number of the L-system (N), and the decay parameter (d). Here, we chose w bc as the changed parameters; the structure will be generated under the defined rules. Because all the input parameters are continuous, the genetic algorithm (GA) is able to 'learn' how to work with the growth behavior and tune it to create better performing designs over time.

Optimization Model and Option Settings
This behavioral generative geometry model can create numerous layout patterns for a wing surface structure based on a relatively small set of input parameters (87 boundary seeds). To find high-performing layout patterns, a series of judgments are needed to the algorithm. In this paper, static finite element analysis (FEA) is used to simulate the performance of each design under the given loading conditions. This analysis gives us a set of metrics that we can use to establish the objectives and constraints of our optimization problem:

•
The total weight of the whole model must be minimum, with the first buckling loads confined to a certain range. • The first buckling loads, which are the objective of this design, should be maximum. These two optimization problems are presented as follows.
where w i is the weight data for the boundary seeds, and M bic is the total weight of the bionic plane. λ RBL is the relative buckling load factor related to the initial straight plane weight (λ RBL = λ·(1 − ((M bic − M initial ))/M initial )). Using the above optimization model, we performed an optimization using a variant of the multi-island genetic algorithm with the following settings in Table 2.  Figure 14 shows the improvement process under different objects and load conditions, evidencing a remarkable ability to develop subsequent generations of designs until the optimal design. Figure 15 illustrates the GA points under uniaxial and biaxial compression.

Result Analysis
where is the weight data for the boundary seeds, and M is the total weight of the bionic plane. is the relative buckling load factor related to the initial straight plane weight ( = λ • (1 − ((M − M ))/M )). Using the above optimization model, we performed an optimization using a variant of the multi-island genetic algorithm with the following settings in Table 2.  Figure 14 shows the improvement process under different objects and load conditions, evidencing a remarkable ability to develop subsequent generations of designs until the optimal design. Figure 15 illustrates the GA points under uniaxial and biaxial compression.   In the uniaxial compression condition, two optimized objects, mass and relative buckling load, corresponding to Equations (6) and (7), are applied to obtain the desired result. The final layout is shown in Figure 16a,b. Figure 16c shows the optimized stiffener layout under biaxial compression with a relative buckling load object corresponding to Equation (7).

Result Analysis
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 17 In the uniaxial compression condition, two optimized objects, mass and relative buckling load, corresponding to Equation (6) and (7), are applied to obtain the desired result. The final layout is shown in Figure 16a,b. Figure 16c shows the optimized stiffener layout under biaxial compression with a relative buckling load object corresponding to Equation (7). The results of both the initial straight and bionic optimized stiffeners under different situations are given in Table 3, and their detailed finite element analysis result are given in Figures 17 and 18, respectively. It is clearly seen that the bionic curved stiffeners increase the critical buckling loads greatly (the highest improvements can reach 40.15%), and the different buckling modes reflect that the bionic method in this paper may divide the global buckling area into pieces, while the local buckling damage plays a dominant role in the buckling mode. As local buckling does not indicate a completely broken situation, this stiffener layout also enlarges the damage tolerance to the solar panel.  The results of both the initial straight and bionic optimized stiffeners under different situations are given in Table 3, and their detailed finite element analysis result are given in Figures 17 and 18, respectively. It is clearly seen that the bionic curved stiffeners increase the critical buckling loads greatly (the highest improvements can reach 40.15%), and the different buckling modes reflect that the bionic method in this paper may divide the global buckling area into pieces, while the local buckling damage plays a dominant role in the buckling mode. As local buckling does not indicate a completely broken situation, this stiffener layout also enlarges the damage tolerance to the solar panel. Table 3. Relative buckling loads (RBLs) and mass data of the initial straight and bionic optimal curved stiffeners under uniaxial compression and biaxial comparison. respectively. It is clearly seen that the bionic curved stiffeners increase the critical buckling loads greatly (the highest improvements can reach 40.15%), and the different buckling modes reflect that the bionic method in this paper may divide the global buckling area into pieces, while the local buckling damage plays a dominant role in the buckling mode. As local buckling does not indicate a completely broken situation, this stiffener layout also enlarges the damage tolerance to the solar panel.

Wing Surface Structure Application
In this section, the proposed two-step bionic curved stiffener layout generation method is implemented for the wing surface structural design. As illustrated in Figure 19, a simplified model with boundary and loading conditions is presented. The optimization problem is listed in Equation (8), and the optimization options are shown in Table 4. With the GS data, a wing box is set to simulate the bionic grid wing's structural properties.

Wing Surface Structure Application
In this section, the proposed two-step bionic curved stiffener layout generation method is implemented for the wing surface structural design. As illustrated in Figure 19, a simplified model with boundary and loading conditions is presented. The optimization problem is listed in Equation (8), and the optimization options are shown in Table 4. With the GS data, a wing box is set to simulate the bionic grid wing's structural properties.

Wing Surface Structure Application
In this section, the proposed two-step bionic curved stiffener layout generation method is implemented for the wing surface structural design. As illustrated in Figure 19, a simplified model with boundary and loading conditions is presented. The optimization problem is listed in Equation (8), and the optimization options are shown in Table 4. With the GS data, a wing box is set to simulate the bionic grid wing's structural properties.   Here, θ twist is the wing's maximum twist angle, U tip is the displacement of the wing tip, and S max is the maximum stress value of the stiffeners. The optimized process is presented in Figure 20. In this figure, the color represents the tip displacement that corresponds to the bending stiffness; blue indicates a low level and red indicates a high level. The detailed FEM result are shown in Figure 21. For the upper and bottom surfaces' optimal results in Figure 21b,c, the difference is caused by the opposite load condition to the bottom surface suffered by the tension force, while the upper surface resists the compression force. The layout of the upper surface revealed a large set of stiffeners towards the wingspan, and the bottom surface revealed a chord-wise like layout, which is a twist that resists performance to improve torsional rigidity.
For the upper and bottom surfaces' optimal results in Figure 21b,c, the difference is caused by the opposite load condition to the bottom surface suffered by the tension force, while the upper surface resists the compression force. The layout of the upper surface revealed a large set of stiffeners towards the wingspan, and the bottom surface revealed a chord-wise like layout, which is a twist that resists performance to improve torsional rigidity. A comparison between the bionic structure and traditional structure is listed in Table 5. The three types are under the same load conditions, and the grid structure uses the same cross section. The bionic grid structure has a 29.1% weight loss compared with the traditional structure and a 22.8% weight loss compared with the initial straight structure. In the table, this surface grid structure shows great advantages in weight loss. As the traditional structure is a beam and slab structure, this weight loss is mainly relayed to the materials' improvement. These remarkable reductions are caused by the different uses of the structure's layout, and the novel bionic design of the stiffener can further reduce the weight. However, we must discuss the twist angle along the wing span ( Figure 22). This structure has a discontinuous torsional rigidity, which may be related to the continuity verification constraint in the optimization. A comparison between the bionic structure and traditional structure is listed in Table 5. The three types are under the same load conditions, and the grid structure uses the same cross section. The bionic grid structure has a 29.1% weight loss compared with the traditional structure and a 22.8% weight loss compared with the initial straight structure. In the table, this surface grid structure shows great advantages in weight loss. As the traditional structure is a beam and slab structure, this weight loss is mainly relayed to the materials' improvement. These remarkable reductions are caused by the different uses of the structure's layout, and the novel bionic design of the stiffener can further reduce the weight. However, we must discuss the twist angle along the wing span ( Figure 22). This structure has a discontinuous torsional rigidity, which may be related to the continuity verification constraint in the optimization. Table 5. Property comparison between the bionic structure and traditional structure.

Design Object and Constraint
Optimal Result

Graphical representation
A comparison between the bionic structure and traditional structure is listed in Table 5. The three types are under the same load conditions, and the grid structure uses the same cross section. The bionic grid structure has a 29.1% weight loss compared with the traditional structure and a 22.8% weight loss compared with the initial straight structure. In the table, this surface grid structure shows great advantages in weight loss. As the traditional structure is a beam and slab structure, this weight loss is mainly relayed to the materials' improvement. These remarkable reductions are caused by the different uses of the structure's layout, and the novel bionic design of the stiffener can further reduce the weight. However, we must discuss the twist angle along the wing span ( Figure 22). This structure has a discontinuous torsional rigidity, which may be related to the continuity verification constraint in the optimization. A comparison between the bionic structure and traditional structure is listed in Table 5. The three types are under the same load conditions, and the grid structure uses the same cross section. The bionic grid structure has a 29.1% weight loss compared with the traditional structure and a 22.8% weight loss compared with the initial straight structure. In the table, this surface grid structure shows great advantages in weight loss. As the traditional structure is a beam and slab structure, this weight loss is mainly relayed to the materials' improvement. These remarkable reductions are caused by the different uses of the structure's layout, and the novel bionic design of the stiffener can further reduce the weight. However, we must discuss the twist angle along the wing span ( Figure 22). This structure has a discontinuous torsional rigidity, which may be related to the continuity verification constraint in the optimization. A comparison between the bionic structure and traditional structure is listed in Table 5. The three types are under the same load conditions, and the grid structure uses the same cross section. The bionic grid structure has a 29.1% weight loss compared with the traditional structure and a 22.8% weight loss compared with the initial straight structure. In the table, this surface grid structure shows great advantages in weight loss. As the traditional structure is a beam and slab structure, this weight loss is mainly relayed to the materials' improvement. These remarkable reductions are caused by the different uses of the structure's layout, and the novel bionic design of the stiffener can further reduce the weight. However, we must discuss the twist angle along the wing span ( Figure 22). This structure has a discontinuous torsional rigidity, which may be related to the continuity verification constraint in the optimization.  Figure 22. The curve of the torsion angle along the wing-span for the three kinds of structures.

Conclusions
In this paper, a new concept of bionic curved grid-stiffened composite structures is proposed for a wing grid-shell composite structure design. Drawing on the biological structures of nature, a stepped automated method was shown in the paper to simulate the growth mode of a stiffener for the whole design domain. In this way, desired sparsely distributed stiffeners with more flexible skin are proposed; these stiffeners are also more effective than traditional ones. The proposed methodology also established a design system for a plane-like structure that can be applied to a beamlike internal structure design, taking advantage of the formal freedom facilitated by recent advances in additive manufacturing. Finally, this nature-based generative design system offers the possibility to be combined with structural data and the bionic generative method to leverage the power of evolutionary computing to derive unique, high-performing solutions for complex design challenges.

Conflicts of Interest:
The authors declare no conflict of interest. Figure 22. The curve of the torsion angle along the wing-span for the three kinds of structures.

Conclusions
In this paper, a new concept of bionic curved grid-stiffened composite structures is proposed for a wing grid-shell composite structure design. Drawing on the biological structures of nature, a stepped automated method was shown in the paper to simulate the growth mode of a stiffener for the whole design domain. In this way, desired sparsely distributed stiffeners with more flexible skin are proposed; these stiffeners are also more effective than traditional ones. The proposed methodology also established a design system for a plane-like structure that can be applied to a beam-like internal structure design, taking advantage of the formal freedom facilitated by recent advances in additive manufacturing. Finally, this nature-based generative design system offers the possibility to be combined with structural data and the bionic generative method to leverage the power of evolutionary computing to derive unique, high-performing solutions for complex design challenges.