# Multidiscipline Topology Optimization of Stiffened Plate/Shell Structures Inspired by Growth Mechanisms of Leaf Veins in Nature

Authors: Baotong Li; Jun Hong; Suna Yan; Zhifeng Liu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Biological structures with preeminent performance in nature endow inexhaustible inspiration for creative design in engineering. In this paper, based on the observation of the natural morphogenesis of leaf veins, we put forward a simple and practical multidiscipline topology optimization method to produce the stiffener layout for plate/shell structures. This method simulates the emergence of complex branching patterns copying the self-optimization of leaf veins which always try to grow into a configuration with global optimal performances. Unlike the conventional topology optimization methods characterized by “subtraction mode,” the proposed method is based on the “addition mode,” giving great potential for designers to achieve more clear stiffener layout patterns rather than vague material distributions and, consequently, saving computational resources as well as enhancing availability of design outputs. Numerical studies of both static and dynamic problems considered in this paper clearly show the suitability of the proposed method for the optimal design of stiffened plate/shell structures.


---

## Body

## 1. Introduction

To satisfy the ever-increasing demands of resource and energy savings, products with light-weight design are desirable in manufacturing industries. One common and cost-effective approach is the application of stiffened plate/shell structures. For example, load-bearing components in large-scale equipments, such as the bed of machine tool, crossbeam of crane, ship hull, and automobile body, are usually box structures filled with stiffener plates distributed horizontally and vertically. It is no doubt that the mechanical performances of the load-bearing components, including stiffness, strength, and dynamic characteristics, are dependent on the layout pattern of internal stiffener plates to a great extent; thus considerable efforts have been devoted to the optimal design of stiffer layout patterns.Following the advancement of the numerical calculation technique, the structural topology optimization design has gotten fast development over the last two decades, offering great potential for better solutions to stiffener layout design. Luo and Gea [1, 2] developed a systematic topology optimization approach to assist the optimal stiffener design of a 3D plate/shell structure, in which both stiffener location and orientation were considered to deal with static as well as eigenvalue problems. Krog and Olhoff [3] dealt with topology optimization of externally rib-stiffened and internally stiffened honey comb and sandwich Mindlin plates by an improved homogenization approach. Bojczuk and Szteleblak [4] presented and validated a computationally efficient method to calculate the sensitivity of stiffener introduction with respect to arbitrary objective functional of displacements, strains, stresses, and resonant frequencies, which can be used to determine the optimum stiffener location. Ansola et al. [5] proposed a global optimality criteria method for the presentation of both the geometry of the shell mid-plane and the layout of surface stiffeners on the shell structure. Other approaches and contributions can be found in the literature [6–10].Since the time when the idea of stiffener layout design has been transformed into a searching problem of optimal material distribution, the subject has been continuously developed for more and more complicated structures, such as the packaging structures [11], hull structures [12, 13], and aeronautical structures [14, 15]. The previous engineering practices in topology optimization can be deemed as a “subtraction-mode-” based procedure, because the resulting layout is extracted from the design domain by eliminating redundant material iteratively. However, such kind of subtraction mode approaches can only produce a vague material distribution; the exact information about the stiffener location and orientation is not available. Therefore, expensive computational efforts and additional postprocessing treatment are required which make the design process tedious, not effective, and time consuming.To overcome these difficulties, it would be a meaningful way to borrow some experiences from nature. Using bio-inspired approach to improve the quality of structure design has attracted increasing attention in academic and engineering field. Actually, problems on how to define the optimal load-bearing topologies are not confined to those mechanical components in machine tool, airplane, or automobile. Elaborate structures in organs in living system also provide interesting problems such as the growth of leaf veins which enable leaves to withstand the self-weight and environmental loads. It would be very natural to presume that biological configurations are formed almost deterministically so that load-bearing functions needed in organs are created efficiently. With the inspiration of branching patterns in nature, the optimal design of stiffener layout can be interpreted as an analog to a growth process of leaf veins, in which the genes of plant leaf yield a potential capability for their tissue cells to branch and to degenerate adaptively so as to grow into a steady self-optimum structure. In this end, the emergence of stiffener layout for engineering design is definitely an “addition-mode-” based procedure that is completely opposite to the subtraction mode employed in conventional approaches.Stiffener optimal design based on the imitation of branching patterns in nature has first been introduced by Ding and Yamazaki in 2004 [16], and since then, is extensively applied [17–19]. The main idea is inspired by the observation that natural branching patterns such as the leaf venation, root, and axis system of plants always grow in such a manner that global optimal performances can be achieved. Although Ding’s method can be easily implemented, its application depends primarily on an empirical formula, in which several key parameters are selected based on designer’s experience. This is not very helpful especially when searching for the best design in the case of multiobjective and multiconstraint scenarios. Questions associated with the optimization parameters have been left basically not answered. So there is a clear need to develop models that can expand and deepen our understanding of the principles, properties, and mechanisms of the emergence of branching patterns in natural systems and that can also foster simple and practical methodologies for engineering design. In this regard, the present study explores the analytical laws and their corresponding equations that underlie morphogenesis phenomena in nature so as to devise a paradigm to simulate the generation process of branching patterns for topology optimization in engineering.The paper is organized as follows. After this introduction, we briefly present the elements of branching patterns of leaf veins and their geometrical characteristics which are pertinent to the present study. Following this description, a well-founded mathematical explanation for the adaptive growth of branching patterns is derived from the Kuhn-Tucker conditions, leading to a novel optimality criterion that can serve engineering purpose for stiffener layout design. With the suggested method, we illustrate the application in static, dynamic, and multidiscipline design problems. Finally, we conclude the paper with the main findings of the present work.


## 2. Parametric Characterization of Leaf Veins in Nature

### 2.1. Geometrical Modeling of Leaf Veins

After billions of years of evolution, nature has created complex and ingenious topologies for solving engineering problems and improving design methods. A remarkable class of such topologies is the branching configuration existed in numerous biological systems whose evolution and function have fascinated and engaged scientists and engineers for centuries. Despite the advancement of computational techniques, the process of finding analytical laws that represent the growth process of branching patterns in nature still remains uncertain. Reproduction of biological growth faces challenging problems as the experimental data can only tell us what is happening, but not why it happens, and the databases of simulation results tend to be dynamic, incomplete, redundant, and very large.Since the general approach to modeling such kind of common growth mechanism is not expected at the present stage, it would make sense to conduct case studies for some typical prototypes, in which their geometrical configurations and physical functions are described mathematically and certain relations are discussed between them. The leaf venation treated here is a living organ that not only possesses the physiological functions, but also adapts to the complex environmental loads by evolving itself into an optimum structure with two major structural patterns: (1) medial axis pattern and (2) closed loop pattern, as shown in Figure1.Figure 1
Leaf vein’s branching patterns in different scales.The morphogenesis of leaf veins can be deemed as (1) a continuous growth process including sprouting, branching, and degenerating activities which can finally form a complicated dichotomous hierarchical system; (2) a continuous optimization process in which both the growth rate and growth direction are adjusted dynamically so as to make it possible for the layout of leaf venation to the dependent on the growth environment. In order to apply the growth optimality to practical engineering design, our first priority is to establish a refined leaf model with venation systems so that the biological growth procedure can be simulated in a more realistic way.Although the morphologies and the growth environments of various leaf veins differ from each other, the geometrical characteristics of leaf venation can be commonly classified into two categories: one is the geometrical properties of the leaf vein itself, that is, its location, orientation, length, and thickness; the other is those of a region in the leaf blade governed by a branch, that is, its contour shape and volume. Based on this, the growth of leaf venation can be simulated as a generation process for stiffeners on a base plate as shown in Figure2.Figure 2
Geometrical modeling of leaf venation.It can be found that the geometrical model includes both the base plate and candidate stiffeners. The base plate is utilized to model a certain region in the leaf blade governed by a new branch, the material modulus of which is defined asEb. The plate center is assumed as the sprouting point, around which several candidate stiffeners are located symmetrically representing various possible sprouting directions of the new branch. The material modulus of the candidate stiffener is denoted as Es. In this preliminary study, we first fix the value of Es/Eb at 1, and the cross-sectional shape of the stiffener is supposed to be a rectangle whose width (ts) is the same as the thickness of the base plate (hb). Therefore, the weight of the stiffener, which will be regarded as a design variable in the simulation process, is changed only depending on the height of the stiffener (hs).


### 2.2. Numerical Observation of Growth Competition

As is known, competition among candidate stiffeners during growth takes place by nonoptimal stiffeners being defeated and finally eliminated. Before the calculation continues, several assumptions are made with regard to simulating the growth competition, and they are as follows (1) an arbitrary branching point along the leaf vein is selected as a sprouting point, around which the local leaf blade is modeled as a simply supported square plate with a ratio of length to thickness oflb/hb=1000. (2) The environmental loads are illustrated in Figure 3, where the base plate is divided into four regions, and two diagonal regions are applied with the same magnitude of the pressure. The ratio of the lower pressure to the higher pressure is defined as 0.2. (3) A certain kind of nutrition is presumed to be distributed over the base plate, which can be deemed as the distribution of the auxin carrier proteins in the mesophyll when a new leaf vein began to grow. In this model, the spatial distribution of the strain energy is simulated as the nutrient distribution, driving the new vein to sprout and branch. (4) A very small initial height (hs0=0.01 mm) is introduced so as to create the seedling stiffeners, whose stiffness is a very small fraction of the whole stiffness and confirmed to have little effect on the structural performance of the base plate in the initial growth period.Figure 3
Growth competition among candidate stiffeners.It should be noted that those parameters defined above may not be accurate, yet the purpose of these assumptions is to reveal the secret design rules derived in nature. In this section, a parametric study focusing on the relationship between the height of candidate stiffeners and the plate strain energy is performed by numerical calculation based on the finite element method.In Figure3, the palettes of changing histories of the plate strain energy are shown as different candidate stiffeners grow from little to mature. It can be found that with the increase in the height of candidate stiffeners, the strain energy of the base plate appears as a decreasing tendency which means that stiffeners grow taller by absorbing the nutrients distributed surrounding. Then, with a further increase in the stiffener height, the average strain energy will eventually reach a convergence which implies that stiffeners are not able to grow because of lack of nutrition. That is, the growth process will stop when the strain energy of the base plate cannot be decreased anymore. Table 1 lists the quantitative relationship between hs and R for each candidate stiffener, where hs refers to the convergence value of the stiffener height, and R refers to the ratio of strain energy between the initial and convergence states.Table 1
Convergence values ofhs and R for each stiffener.
B-1
B-2
B-3
B-4
B-5
B-6
B-7
B-8
h
s (mm)
7.460
9.700
7.460
9.920
7.460
9.700
7.460
9.920
R (%)
31.74
32.03
31.74
33.16
31.74
32.03
31.74
33.16From this table we can approximate where the minimum strain energy will occur as the height of the stiffener increases with the step size. An interesting phenomenon is that although the increase in each stiffener height can eventually lead to a convergence of the strain energy, these stiffeners differ from each other in the convergence values. For instance, when the height of the stiffener_1 increases to 7.46 mm, the strain energy ratio will converge to 31.74%. However, the convergence values of the stiffener height and strain energy ratio are 9.70, 32.03% and 9.92, 33.16% for the stiffener_2 and the stiffener_4, respectively. Such kind of difference reflects the growth competition among those candidate stiffeners. A reasonable view is that there may exist an ideal “balanced point,” at which the weights of the stiffeners follow a distribution in such a way that the average strain energy will be minimized. The more materials the stiffeners obtain, the more opportunities they will get to win in the fierce growth competition, and the final winner will represent the sprouting direction of the new vein.


## 2.1. Geometrical Modeling of Leaf Veins

After billions of years of evolution, nature has created complex and ingenious topologies for solving engineering problems and improving design methods. A remarkable class of such topologies is the branching configuration existed in numerous biological systems whose evolution and function have fascinated and engaged scientists and engineers for centuries. Despite the advancement of computational techniques, the process of finding analytical laws that represent the growth process of branching patterns in nature still remains uncertain. Reproduction of biological growth faces challenging problems as the experimental data can only tell us what is happening, but not why it happens, and the databases of simulation results tend to be dynamic, incomplete, redundant, and very large.Since the general approach to modeling such kind of common growth mechanism is not expected at the present stage, it would make sense to conduct case studies for some typical prototypes, in which their geometrical configurations and physical functions are described mathematically and certain relations are discussed between them. The leaf venation treated here is a living organ that not only possesses the physiological functions, but also adapts to the complex environmental loads by evolving itself into an optimum structure with two major structural patterns: (1) medial axis pattern and (2) closed loop pattern, as shown in Figure1.Figure 1
Leaf vein’s branching patterns in different scales.The morphogenesis of leaf veins can be deemed as (1) a continuous growth process including sprouting, branching, and degenerating activities which can finally form a complicated dichotomous hierarchical system; (2) a continuous optimization process in which both the growth rate and growth direction are adjusted dynamically so as to make it possible for the layout of leaf venation to the dependent on the growth environment. In order to apply the growth optimality to practical engineering design, our first priority is to establish a refined leaf model with venation systems so that the biological growth procedure can be simulated in a more realistic way.Although the morphologies and the growth environments of various leaf veins differ from each other, the geometrical characteristics of leaf venation can be commonly classified into two categories: one is the geometrical properties of the leaf vein itself, that is, its location, orientation, length, and thickness; the other is those of a region in the leaf blade governed by a branch, that is, its contour shape and volume. Based on this, the growth of leaf venation can be simulated as a generation process for stiffeners on a base plate as shown in Figure2.Figure 2
Geometrical modeling of leaf venation.It can be found that the geometrical model includes both the base plate and candidate stiffeners. The base plate is utilized to model a certain region in the leaf blade governed by a new branch, the material modulus of which is defined asEb. The plate center is assumed as the sprouting point, around which several candidate stiffeners are located symmetrically representing various possible sprouting directions of the new branch. The material modulus of the candidate stiffener is denoted as Es. In this preliminary study, we first fix the value of Es/Eb at 1, and the cross-sectional shape of the stiffener is supposed to be a rectangle whose width (ts) is the same as the thickness of the base plate (hb). Therefore, the weight of the stiffener, which will be regarded as a design variable in the simulation process, is changed only depending on the height of the stiffener (hs).


## 2.2. Numerical Observation of Growth Competition

As is known, competition among candidate stiffeners during growth takes place by nonoptimal stiffeners being defeated and finally eliminated. Before the calculation continues, several assumptions are made with regard to simulating the growth competition, and they are as follows (1) an arbitrary branching point along the leaf vein is selected as a sprouting point, around which the local leaf blade is modeled as a simply supported square plate with a ratio of length to thickness oflb/hb=1000. (2) The environmental loads are illustrated in Figure 3, where the base plate is divided into four regions, and two diagonal regions are applied with the same magnitude of the pressure. The ratio of the lower pressure to the higher pressure is defined as 0.2. (3) A certain kind of nutrition is presumed to be distributed over the base plate, which can be deemed as the distribution of the auxin carrier proteins in the mesophyll when a new leaf vein began to grow. In this model, the spatial distribution of the strain energy is simulated as the nutrient distribution, driving the new vein to sprout and branch. (4) A very small initial height (hs0=0.01 mm) is introduced so as to create the seedling stiffeners, whose stiffness is a very small fraction of the whole stiffness and confirmed to have little effect on the structural performance of the base plate in the initial growth period.Figure 3
Growth competition among candidate stiffeners.It should be noted that those parameters defined above may not be accurate, yet the purpose of these assumptions is to reveal the secret design rules derived in nature. In this section, a parametric study focusing on the relationship between the height of candidate stiffeners and the plate strain energy is performed by numerical calculation based on the finite element method.In Figure3, the palettes of changing histories of the plate strain energy are shown as different candidate stiffeners grow from little to mature. It can be found that with the increase in the height of candidate stiffeners, the strain energy of the base plate appears as a decreasing tendency which means that stiffeners grow taller by absorbing the nutrients distributed surrounding. Then, with a further increase in the stiffener height, the average strain energy will eventually reach a convergence which implies that stiffeners are not able to grow because of lack of nutrition. That is, the growth process will stop when the strain energy of the base plate cannot be decreased anymore. Table 1 lists the quantitative relationship between hs and R for each candidate stiffener, where hs refers to the convergence value of the stiffener height, and R refers to the ratio of strain energy between the initial and convergence states.Table 1
Convergence values ofhs and R for each stiffener.
B-1
B-2
B-3
B-4
B-5
B-6
B-7
B-8
h
s (mm)
7.460
9.700
7.460
9.920
7.460
9.700
7.460
9.920
R (%)
31.74
32.03
31.74
33.16
31.74
32.03
31.74
33.16From this table we can approximate where the minimum strain energy will occur as the height of the stiffener increases with the step size. An interesting phenomenon is that although the increase in each stiffener height can eventually lead to a convergence of the strain energy, these stiffeners differ from each other in the convergence values. For instance, when the height of the stiffener_1 increases to 7.46 mm, the strain energy ratio will converge to 31.74%. However, the convergence values of the stiffener height and strain energy ratio are 9.70, 32.03% and 9.92, 33.16% for the stiffener_2 and the stiffener_4, respectively. Such kind of difference reflects the growth competition among those candidate stiffeners. A reasonable view is that there may exist an ideal “balanced point,” at which the weights of the stiffeners follow a distribution in such a way that the average strain energy will be minimized. The more materials the stiffeners obtain, the more opportunities they will get to win in the fierce growth competition, and the final winner will represent the sprouting direction of the new vein.


## 3. Mathematical Insights into Growth Competition

In order to capture some of the most essential features of the growth competition in leaf venation, an optimization concept termed guide-weight criterion [20, 21] is introduced in this section, by means of which the relationship between the growth objective and the corresponding resource constraint will be clarified to serve engineering purpose in stiffened plate/shell structure design.


### 3.1. Modeling of the “Balanced Point”

Actually, there are two ways for a candidate stiffener to survive the competition for best. One is to reduce the amount of material used while maintaining the same mechanical performance; the other is to maintain the amount of the material used but with improved mechanical performance. To transform the implicit “balanced point” into a group of explicit mathematical equations, the Kuhn-Tucker conditions are employed to describe the growth competition problem, shown as follows:(1)findW=[W1,W2,…,WN]T∈RN,min:F(W)=∑j=1Jrj·fj(W),subject  toG(W)≤Wg,WiU>Wi>WiL,i=1,2,…,N,
where W=[W1,W2,…,WN]T is an N-dimensional design vector, and Wi is the weight of the ith candidate stiffener; F(W) is the weighted growth objective function, fj(W) is the jth subobjective function which denotes the mechanical performance such as the strain energy, strength, and natural frequency, rj is the corresponding weight factor; G(W) is the resource constraint function which denotes the amount of the material allocated to candidate stiffeners, and Wg is the upper limit of resource constraint.The Lagrange function is constructed as follows:(2)L(W)=F(W)+λ·(G(W)-Wg),
where λ is the Lagrange multiplier.Based on the Kuhn-Tucker conditions, at the optimal solutionW*, the following formulas must be satisfied:(C-1)(3)∇L=∇F(W)+λ·∇[G(W)-Wg]==={≤0,Wi=WiU,=0,WiU>Wi>WiL,≥0,Wi=WiL,(C-2)(4)λ≥0,(C-3)(5)===λ·[G(W)-Wg]=0,{λ=0,G(W)-Wg<0,λ≥0,G(W)-Wg=0,(C-4)(6)G(W)≤Wg,(C-5)(7)WiU>Wi>WiL,i=1,2,…,N.MultiplyingWi by (3) leads to
(8)Wi·∂F∂Wi+λ·Wi·∂G∂Wi=0,because(9)∂G∂Wi=∂(∑i=1NWi)∂Wi=1.When applying (9), (8) can be rewritten as
(10)Wi·∂F∂Wi+λ·Wi=0.Hence, the weight of theith candidate stiffener can be expressed as
(11)Wi=-1λ·(Wi·∂F∂Wi)=-1λ·[Wi·(∑j=1Jrj·∂fj∂Wi)].The total weight of the stiffeners can be gained from the following equation:(12)Wsum=∑i=1NWi=-1λ·∑i=1N[Wi·(∑j=1Jrj·∂fj∂Wi)].Let(13)Di=-Wi·(∑j=1Jrj·∂fj∂Wi),(14)Dsum=∑i=1N[-Wi·(∑j=1Jrj·∂fj∂Wi)],
where Di is the generalized weight of the ith candidate stiffener which is the negative product of the stiffener weight and the sensitivity of the growth objective function with respect to the stiffener weight, and Dsum is the generalized total weight of all the candidate stiffeners.Thus, the weight of theith candidate stiffener can be expressed by a simplified form:
(15)Wi=(Di·WsumDsum),(i=1,…,N).When iterating, it can be written as(16)Wi(k+1)=(Di·WsumDsum)(k),(i=1,…,N).The above equation is the so-called guide-weight criterion that yields an opportunity for designers to get the ideal “balanced point” in terms of weight distribution of candidate stiffeners.In order to guarantee the convergence of the calculation, a relaxation factorα is introduced, and the iterative forms can be rewritten as
(17)Wi(k+1)=α(k)·(Di·WsumDsum)(k)+(1-α(k))·Wi(k),------------------(i=1,…,N).From (17), we establish the governing criterion to allocate the growth resources to candidate stiffeners, guiding the growth process to make the weight distribution satisfy the previous Kuhn-Tucker conditions. Based on the observation and derivations above, the “balanced point” among candidate stiffeners can be assumed to take the following concise form:
(18)D1W1=D2W2=⋯DiWi=⋯=DsumWsum=λ,==============-=(i=1,…,N).Equation (18) shows that in an optimal growth result, the material allocated to candidate stiffeners should be exactly proportional to the value of their generalized weight, and the “balanced point” is the Lagrange multiplier λ. That is, on the premise of unchanged growth resources, the growth result is optimal when the positive change of the objective function caused by the weight increase of any one of the involved candidate stiffeners is identical to the negative change caused by the weight decrease of any other candidate stiffeners.The above-mentioned “balanced point” can be evaluated quantitatively by the following equation:(19)σ=∑i=1N(Wi/WsumDi/Dsum-1)2,(i=1,…,N),
where σ is called the “balanced degree” which should be zero for an optimal structure.In this paper, the “balanced degree” is utilized to determine the step factorα(k):
(20)α(k)=lne3N·ln(σ(k)-0.1σ(0)+1),
where σ(0) is the initial value of the “balanced degree,” σ(k) is the “balanced degree” after kth iteration, and N is the total number of the candidate stiffeners. The iteration process would be terminated when the value of σ(k) is decreased to be one-tenth of its initial value σ(0).


### 3.2. Mathematical Explanation of the “Balanced Point”

The findings in the previous section motivate us to explore an explanation, where the mathematical nature about the relationship among the Lagrange multiplierλ, objective function F(X), and constraint boundary Wg will be clarified, so that the natural growth mechanism could be fully employed to engineering design.By taking the first derivative of (3) with respect to the constraint boundary Wg, we have the following relation:
(21)∂L∂Wg=∂F∂Wg+∂λ∂Wg(G-Wg)+λ·∂G∂Wg-λ.Let(22)∇WgW=[∂W1∂Wg,∂W2∂Wg,…,∂WN∂Wg]T.Equation (21) can be rewritten as
(23)∂L∂Wg=∇WgW·(∂F∂W+λ·∂G∂W)=+∂λ∂Wg(G-Wg)-λ,since:(24)∂F∂W+λ·∂G∂W=0.Equation (21) can be further simplified as
(25)∂L∂Wg=∂λ∂Wg(G-Wg)-λ.According to (5), we obtain
(26)∂λ∂Wg(G-Wg)=0.Thus(27)∂L∂Wg=-λ.Consider that the optimal solutionW* satisfies the following equation:
(28)L(W*)=F(W*).Thus, the relationship among the Lagrange multiplierλ, objective function F(X), and constraint boundary Wg can be expressed by
(29)∂F∂Wg=-λ.The above equation provides very important insights into the essential features of the growth mechanism which is most efficient in carrying biological structures, such as leaf veins. Note that the negative Lagrange multiplier is equal to the sensitivity ofF(X) with respect to Wg; hence there are three cases, that is, λ>0, λ=0, and λ<0, during the iterative growth process. (1)
Whenλ>0, according to (18), one can find that Di > 0, that is, ∂F/∂Wi<0, which means the weight of the ith candidate stiffener Wi should be increased with larger Di so as to get an optimal growth result. That is, the objective function F(X) would be improved with an increase in constraint boundary Wg.(2)
Whenλ<0, it can be seen that the constraint boundary Wg should be decreased in order to improve the objective function F(X). That means that the total weight of all the candidate stiffeners Wsum should be decreased with larger λ after each iterative growth step until λ=0, so that the Kuhn-Tucker conditions would be satisfied.(3)
Whenλ=0, whether for an increase or for a decrease in the constraint boundary Wg, the objective function F(X) cannot be improved anymore.


## 3.1. Modeling of the “Balanced Point”

Actually, there are two ways for a candidate stiffener to survive the competition for best. One is to reduce the amount of material used while maintaining the same mechanical performance; the other is to maintain the amount of the material used but with improved mechanical performance. To transform the implicit “balanced point” into a group of explicit mathematical equations, the Kuhn-Tucker conditions are employed to describe the growth competition problem, shown as follows:(1)findW=[W1,W2,…,WN]T∈RN,min:F(W)=∑j=1Jrj·fj(W),subject  toG(W)≤Wg,WiU>Wi>WiL,i=1,2,…,N,
where W=[W1,W2,…,WN]T is an N-dimensional design vector, and Wi is the weight of the ith candidate stiffener; F(W) is the weighted growth objective function, fj(W) is the jth subobjective function which denotes the mechanical performance such as the strain energy, strength, and natural frequency, rj is the corresponding weight factor; G(W) is the resource constraint function which denotes the amount of the material allocated to candidate stiffeners, and Wg is the upper limit of resource constraint.The Lagrange function is constructed as follows:(2)L(W)=F(W)+λ·(G(W)-Wg),
where λ is the Lagrange multiplier.Based on the Kuhn-Tucker conditions, at the optimal solutionW*, the following formulas must be satisfied:(C-1)(3)∇L=∇F(W)+λ·∇[G(W)-Wg]==={≤0,Wi=WiU,=0,WiU>Wi>WiL,≥0,Wi=WiL,(C-2)(4)λ≥0,(C-3)(5)===λ·[G(W)-Wg]=0,{λ=0,G(W)-Wg<0,λ≥0,G(W)-Wg=0,(C-4)(6)G(W)≤Wg,(C-5)(7)WiU>Wi>WiL,i=1,2,…,N.MultiplyingWi by (3) leads to
(8)Wi·∂F∂Wi+λ·Wi·∂G∂Wi=0,because(9)∂G∂Wi=∂(∑i=1NWi)∂Wi=1.When applying (9), (8) can be rewritten as
(10)Wi·∂F∂Wi+λ·Wi=0.Hence, the weight of theith candidate stiffener can be expressed as
(11)Wi=-1λ·(Wi·∂F∂Wi)=-1λ·[Wi·(∑j=1Jrj·∂fj∂Wi)].The total weight of the stiffeners can be gained from the following equation:(12)Wsum=∑i=1NWi=-1λ·∑i=1N[Wi·(∑j=1Jrj·∂fj∂Wi)].Let(13)Di=-Wi·(∑j=1Jrj·∂fj∂Wi),(14)Dsum=∑i=1N[-Wi·(∑j=1Jrj·∂fj∂Wi)],
where Di is the generalized weight of the ith candidate stiffener which is the negative product of the stiffener weight and the sensitivity of the growth objective function with respect to the stiffener weight, and Dsum is the generalized total weight of all the candidate stiffeners.Thus, the weight of theith candidate stiffener can be expressed by a simplified form:
(15)Wi=(Di·WsumDsum),(i=1,…,N).When iterating, it can be written as(16)Wi(k+1)=(Di·WsumDsum)(k),(i=1,…,N).The above equation is the so-called guide-weight criterion that yields an opportunity for designers to get the ideal “balanced point” in terms of weight distribution of candidate stiffeners.In order to guarantee the convergence of the calculation, a relaxation factorα is introduced, and the iterative forms can be rewritten as
(17)Wi(k+1)=α(k)·(Di·WsumDsum)(k)+(1-α(k))·Wi(k),------------------(i=1,…,N).From (17), we establish the governing criterion to allocate the growth resources to candidate stiffeners, guiding the growth process to make the weight distribution satisfy the previous Kuhn-Tucker conditions. Based on the observation and derivations above, the “balanced point” among candidate stiffeners can be assumed to take the following concise form:
(18)D1W1=D2W2=⋯DiWi=⋯=DsumWsum=λ,==============-=(i=1,…,N).Equation (18) shows that in an optimal growth result, the material allocated to candidate stiffeners should be exactly proportional to the value of their generalized weight, and the “balanced point” is the Lagrange multiplier λ. That is, on the premise of unchanged growth resources, the growth result is optimal when the positive change of the objective function caused by the weight increase of any one of the involved candidate stiffeners is identical to the negative change caused by the weight decrease of any other candidate stiffeners.The above-mentioned “balanced point” can be evaluated quantitatively by the following equation:(19)σ=∑i=1N(Wi/WsumDi/Dsum-1)2,(i=1,…,N),
where σ is called the “balanced degree” which should be zero for an optimal structure.In this paper, the “balanced degree” is utilized to determine the step factorα(k):
(20)α(k)=lne3N·ln(σ(k)-0.1σ(0)+1),
where σ(0) is the initial value of the “balanced degree,” σ(k) is the “balanced degree” after kth iteration, and N is the total number of the candidate stiffeners. The iteration process would be terminated when the value of σ(k) is decreased to be one-tenth of its initial value σ(0).


## 3.2. Mathematical Explanation of the “Balanced Point”

The findings in the previous section motivate us to explore an explanation, where the mathematical nature about the relationship among the Lagrange multiplierλ, objective function F(X), and constraint boundary Wg will be clarified, so that the natural growth mechanism could be fully employed to engineering design.By taking the first derivative of (3) with respect to the constraint boundary Wg, we have the following relation:
(21)∂L∂Wg=∂F∂Wg+∂λ∂Wg(G-Wg)+λ·∂G∂Wg-λ.Let(22)∇WgW=[∂W1∂Wg,∂W2∂Wg,…,∂WN∂Wg]T.Equation (21) can be rewritten as
(23)∂L∂Wg=∇WgW·(∂F∂W+λ·∂G∂W)=+∂λ∂Wg(G-Wg)-λ,since:(24)∂F∂W+λ·∂G∂W=0.Equation (21) can be further simplified as
(25)∂L∂Wg=∂λ∂Wg(G-Wg)-λ.According to (5), we obtain
(26)∂λ∂Wg(G-Wg)=0.Thus(27)∂L∂Wg=-λ.Consider that the optimal solutionW* satisfies the following equation:
(28)L(W*)=F(W*).Thus, the relationship among the Lagrange multiplierλ, objective function F(X), and constraint boundary Wg can be expressed by
(29)∂F∂Wg=-λ.The above equation provides very important insights into the essential features of the growth mechanism which is most efficient in carrying biological structures, such as leaf veins. Note that the negative Lagrange multiplier is equal to the sensitivity ofF(X) with respect to Wg; hence there are three cases, that is, λ>0, λ=0, and λ<0, during the iterative growth process. (1)
Whenλ>0, according to (18), one can find that Di > 0, that is, ∂F/∂Wi<0, which means the weight of the ith candidate stiffener Wi should be increased with larger Di so as to get an optimal growth result. That is, the objective function F(X) would be improved with an increase in constraint boundary Wg.(2)
Whenλ<0, it can be seen that the constraint boundary Wg should be decreased in order to improve the objective function F(X). That means that the total weight of all the candidate stiffeners Wsum should be decreased with larger λ after each iterative growth step until λ=0, so that the Kuhn-Tucker conditions would be satisfied.(3)
Whenλ=0, whether for an increase or for a decrease in the constraint boundary Wg, the objective function F(X) cannot be improved anymore.


## 4. Biologically Inspired Topology Optimization Algorithm

The preceding modeling and computational efforts focus primarily on the derivation of the optimality criterion that governs the growth behavior of branching patterns in leaf veins and have not yet given an evolutionary programming for the optimal design of natural and engineering systems. This section will elaborate how to employ such optimality criterion to reproduce the growth process of branching patterns in leaf veins and how to solve the multidiscipline topology optimization problem for stiffened plate/shell structures.


### 4.1. Simulation Strategy

The objective of this problem is to seek the optimal structural performance by growing a given amount of material as stiffeners. The reproduction of adaptive growth process is treated by means of a comprehensive iterative three-step approach, as shown in Figure4.Figure 4
Simulation strategy for reproducing the adaptive growth process.Step 1 (initialization).
The sprouting points for candidate stiffeners to grow are initialized and are included in a branching point set{B}; stiffeners extended from the sprouting points are gathered in the corresponding stiffener set {C}.Step 2 (competition).
The competition for growth resources among candidate stiffeners is simulated by performing one or more material allocation substeps, in which the growth rate of each stiffener is controlled based on the weight distribution optimized by the iterative criterion derived in Section3 (i.e., (17)). The height of stiffener cross section (hs) is utilized as a decision variable in the optimization process. As a result, some of the stiffeners in the set {C} will become taller and taller, while the other stiffeners will not. Since the material allocation is treated as an inequality constrained optimization problem, the total increment of stiffener weight obtained in each iterative step, ΔWi(i=1,2,…,s), is not greater than the given constraint boundary Wg.Step 3 (reconstruction).
During the growth process, the stiffener layout is continuously reconstructed by obeying the branching and degenerating rules like those for leaf veins. For example, when a stiffener grows up to a height that is greater than the specified upper limithb, a subtle control is activated endowing this stiffener with capabilities for branching operation, in which two ends of the stiffener are added into the set {B}, and new stiffeners around the new sprouting points are also supplemented into the set {C}. Likewise, if the height of candidate stiffeners is less than the specified lower limit hd, the degenerating operation will be activated, where the involved stiffeners are eliminated from the set {C} and their ends are also deleted from the set {B}. By doing this, the best potential direction can be selected for stiffeners to grow in the next iterative step.Both competition and reconstruction loops are repeated alternatively until there are no any material resources available from the specified growth environment. That is, the growth process will stop when the weight increment of the candidate stiffeners reaches the given upper limitWp. The difference of the growth process between the macroscale pattern (i.e., the medial axis pattern) and the microscale pattern (i.e., the closed loop pattern) is mainly reflected in the final step of the simulation process (i.e., the reconstruction step). For the medial axis pattern, the only winner in each growth step is just the tallest stiffener that would be maintained for branching, while others are eliminated. However, for the closed loop pattern, stiffeners whose heights are greater than the specified threshold value would be eligible to participate in the next round of competition, thus yielding opportunity to form a closed loop in the following growth process.


### 4.2. Numerical Application

To validate the effectiveness of the reported bio-inspired methodology, three numerical examples are presented by generating the optimal stiffener layouts for a simply supported square plate with nonuniformly distributed pressure shown in Figure3. The first example is a static design problem pursuing the stiffest stiffener layout so as to minimize its average strain energy; the second example is a dynamic design problem that aims to maximize the natural frequency of its first mode; and the third example is a multidiscipline design problem in which the importance of each objective is determined by putting in appropriate weights (i.e., ωStatic=0.6 and ωDynamic=0.4). Before getting into the actual design, the sprouting points are specified firstly on the four central points of the plate edges. The total increment of stiffener weight, Wr, is specified as 1.5Wp, where Wp denotes the initial weight of the base plate. The upper limit for stiffeners to branch is set as hb=104hd, where hd is the lower limit for stiffeners to degenerate. Based on these assumptions, the branch-like stiffener layouts with the weight fraction varied from 0.3 to 1.5 are constructed. To facilitate comparison, a conventional algorism based on the density method is also applied to search for the optimal material distribution, in which the ratio of retained material β is controlled from 30% to 70%, as shown in Figures 5, 6, 7, and 8.Figure 5
Comparison of stiffener layouts produced by the suggested method and the density method for static design problem.Figure 6
Comparison of stiffener layouts produced by the suggested method and the density method for dynamic design problem.Figure 7
Comparison of stiffener layouts produced by the suggested method and the density method for multidiscipline design problem (strain energy).Figure 8
Comparison of stiffener layouts produced by the suggested method and the density method for multidiscipline design problem (vibration mode).It can be found that the suggested method is able to produce the effective stiffener layouts adaptively to various design objectives. For the static design problem, the strain energy of the final stiffened structure deceases nearly 1927 times compared to the base plate without stiffener; for the dynamic design problem, the first eigenvalue of the final stiffened structure increases nearly 12 times compared to the base plate without stiffener; and for the multidiscipline design problem, the strain energy of the final stiffened structure decreases nearly 325 times with the first eigenvalue increased about 5 times.It can also be found that the simulated branching patterns have both higher and lower height levels, which can be interpreted to be responsible for the load-bearing properties of leaf veins in nature that the higher branches (i.e., the central veins) bear most of the principal stress, and the lower branches (i.e., the lateral veins) are the tracing of the gradient of shear stress. Obviously, the generated stiffeners are primarily concentrated in the regions which are in qualitative consistence with the optimal material distributions resulted from the conventional algorism. However, the complex branching patterns produced by the suggested method give rise to flexibility and distinctiveness of design outputs. The flexibility here means that the bio-inspired method that could generate optimal load-bearing topologies would also possess a higher tolerance to variations of design objectives which is an unintended benefit for survival in a harsh growth environment. Design results produced by the presented method are relatively clear stiffener distribution patterns rather than vague material distributions, which will help to minimize effort, time and cost during the product development.


## 4.1. Simulation Strategy

The objective of this problem is to seek the optimal structural performance by growing a given amount of material as stiffeners. The reproduction of adaptive growth process is treated by means of a comprehensive iterative three-step approach, as shown in Figure4.Figure 4
Simulation strategy for reproducing the adaptive growth process.Step 1 (initialization).
The sprouting points for candidate stiffeners to grow are initialized and are included in a branching point set{B}; stiffeners extended from the sprouting points are gathered in the corresponding stiffener set {C}.Step 2 (competition).
The competition for growth resources among candidate stiffeners is simulated by performing one or more material allocation substeps, in which the growth rate of each stiffener is controlled based on the weight distribution optimized by the iterative criterion derived in Section3 (i.e., (17)). The height of stiffener cross section (hs) is utilized as a decision variable in the optimization process. As a result, some of the stiffeners in the set {C} will become taller and taller, while the other stiffeners will not. Since the material allocation is treated as an inequality constrained optimization problem, the total increment of stiffener weight obtained in each iterative step, ΔWi(i=1,2,…,s), is not greater than the given constraint boundary Wg.Step 3 (reconstruction).
During the growth process, the stiffener layout is continuously reconstructed by obeying the branching and degenerating rules like those for leaf veins. For example, when a stiffener grows up to a height that is greater than the specified upper limithb, a subtle control is activated endowing this stiffener with capabilities for branching operation, in which two ends of the stiffener are added into the set {B}, and new stiffeners around the new sprouting points are also supplemented into the set {C}. Likewise, if the height of candidate stiffeners is less than the specified lower limit hd, the degenerating operation will be activated, where the involved stiffeners are eliminated from the set {C} and their ends are also deleted from the set {B}. By doing this, the best potential direction can be selected for stiffeners to grow in the next iterative step.Both competition and reconstruction loops are repeated alternatively until there are no any material resources available from the specified growth environment. That is, the growth process will stop when the weight increment of the candidate stiffeners reaches the given upper limitWp. The difference of the growth process between the macroscale pattern (i.e., the medial axis pattern) and the microscale pattern (i.e., the closed loop pattern) is mainly reflected in the final step of the simulation process (i.e., the reconstruction step). For the medial axis pattern, the only winner in each growth step is just the tallest stiffener that would be maintained for branching, while others are eliminated. However, for the closed loop pattern, stiffeners whose heights are greater than the specified threshold value would be eligible to participate in the next round of competition, thus yielding opportunity to form a closed loop in the following growth process.


## 4.2. Numerical Application

To validate the effectiveness of the reported bio-inspired methodology, three numerical examples are presented by generating the optimal stiffener layouts for a simply supported square plate with nonuniformly distributed pressure shown in Figure3. The first example is a static design problem pursuing the stiffest stiffener layout so as to minimize its average strain energy; the second example is a dynamic design problem that aims to maximize the natural frequency of its first mode; and the third example is a multidiscipline design problem in which the importance of each objective is determined by putting in appropriate weights (i.e., ωStatic=0.6 and ωDynamic=0.4). Before getting into the actual design, the sprouting points are specified firstly on the four central points of the plate edges. The total increment of stiffener weight, Wr, is specified as 1.5Wp, where Wp denotes the initial weight of the base plate. The upper limit for stiffeners to branch is set as hb=104hd, where hd is the lower limit for stiffeners to degenerate. Based on these assumptions, the branch-like stiffener layouts with the weight fraction varied from 0.3 to 1.5 are constructed. To facilitate comparison, a conventional algorism based on the density method is also applied to search for the optimal material distribution, in which the ratio of retained material β is controlled from 30% to 70%, as shown in Figures 5, 6, 7, and 8.Figure 5
Comparison of stiffener layouts produced by the suggested method and the density method for static design problem.Figure 6
Comparison of stiffener layouts produced by the suggested method and the density method for dynamic design problem.Figure 7
Comparison of stiffener layouts produced by the suggested method and the density method for multidiscipline design problem (strain energy).Figure 8
Comparison of stiffener layouts produced by the suggested method and the density method for multidiscipline design problem (vibration mode).It can be found that the suggested method is able to produce the effective stiffener layouts adaptively to various design objectives. For the static design problem, the strain energy of the final stiffened structure deceases nearly 1927 times compared to the base plate without stiffener; for the dynamic design problem, the first eigenvalue of the final stiffened structure increases nearly 12 times compared to the base plate without stiffener; and for the multidiscipline design problem, the strain energy of the final stiffened structure decreases nearly 325 times with the first eigenvalue increased about 5 times.It can also be found that the simulated branching patterns have both higher and lower height levels, which can be interpreted to be responsible for the load-bearing properties of leaf veins in nature that the higher branches (i.e., the central veins) bear most of the principal stress, and the lower branches (i.e., the lateral veins) are the tracing of the gradient of shear stress. Obviously, the generated stiffeners are primarily concentrated in the regions which are in qualitative consistence with the optimal material distributions resulted from the conventional algorism. However, the complex branching patterns produced by the suggested method give rise to flexibility and distinctiveness of design outputs. The flexibility here means that the bio-inspired method that could generate optimal load-bearing topologies would also possess a higher tolerance to variations of design objectives which is an unintended benefit for survival in a harsh growth environment. Design results produced by the presented method are relatively clear stiffener distribution patterns rather than vague material distributions, which will help to minimize effort, time and cost during the product development.


## 5. Conclusions

In this work, we introduced a bio-inspired topology optimization method that integrates the morphogenesis law of leaf veins with conditioned extremum optimization theory for modeling, generating, and evolving complex branching patterns so as to identify the optimal load-bearing topology of stiffened plate/shell structures. Based on the results presented in this paper, the following conclusions can be drawn.(1)
The reported bio-inspired methodology is a valid topology optimizer by which the design solutions produced are robust and distinct stiffener layout patterns rather than vague material distributions, avoiding the difficulties caused by conventional methods that normally require an additional postprocessing treatment to distinguish the real layout solution. Consequently, design efficiency is substantially improved which will help to minimize effort, time, and cost during the product development.(2)
The present work yields very important insights into the theoretical implications of natural growth mechanism for practical engineering design. Among the key problems currently facing biologists and engineers is how to understand the self-optimization characteristics of branching patterns in nature, such as the geometrical form and physical function causality in the evolution of leaf veins. The proposed simulation strategy can be used to verify postulations in a sufficiently controlled environment, where it is likely to isolate various physical mechanisms of interest.The work presented in this paper is by no means complete. Future studies will include continuous model refinement. For instance, experimental and numerical studies of the relationship between typical distribution patterns of leaf veins and their mechanics self-adaptability will be carried out. The influence of different growth objectives and resource constraints will be thoroughly investigated and correlated with various natural conditions. The simulation of biological growth is a very complex task due to the strong randomness and nonlinear factors associated with the environmental conditions. In order to capture the essential features of biological growth in nature, a robust design optimization model is currently under development so that both the optimality and robustness of the simulated growth results would be taken into account in the future research.



---

*Source: 653895-2013-05-16.xml*