Simultaneous optimization of topology and layout of modular stiffeners on shells and plates

Stiffened shells and plates are widely used in engineering. Their performance is highly influenced by the arrangement, or layout, of stiffeners on the base shell or plate and the geometric features, or topology, of these stiffeners. Moreover, modular design is beneficial, since it allows for increased quality control and mass production. In this work, a method is developed that simultaneously optimizes the topology of stiffeners and their layout on a base shell or plate. This is accomplished by introducing a fixed number of modular stiffeners, which are subject to density-based topology optimization and a mapping of these modules to a ground structure. To illustrate potential applications, several stiffened plates and shell examples are presented. All examples demonstrated that the proposed method is able to generate clear topologies for any number of modules and a distinct layout of the stiffeners on the base shell or plate.


Introduction
Stiffened shells are used widely in, for example, maritime, aerospace, and civil structures. This is because of their high load carrying capacity and light-weight properties. However, due to their thin-walled features, these structures are usually sensitive to loads leading to out-of-plane deflections, imperfections, vibrations, and buckling. Such responses are influenced by the geometric proportions of the stiffeners and base shell, and the layout of the stiffeners on the base shell (Bedair 1997). However, changing the geometric features, such as the thickness of the base shell, as well as every individual stiffener from point to point, is often infeasible due to manufacturing difficulties and costs (Chung and Lee 1997;Lam and Santhikumar 2003). Moreover, every unique component has to be produced and qualified apiece. Thus, the tendency in industry is towards designing structures with fewer component types, since it allows for increased and cheaper quality control, more accessible mass production, and therewith reduction of costs (Higginson et al. 2020). This reuse of components is called modularity, where a module is defined as a component with particular geometric features that can be repeatedly used in the design domain.
The complexity of the above design problem makes the result highly dependent on the experience of the designer. Therefore, a model-based structural optimization technique known as topology optimization poses a solution. Topology optimization has shown the ability to solve complex structural design problems and to produce competitive and innovative solutions (Koppen 2017).
The topology optimization problem that aims to find the optimal layout for stiffeners on a base shell has been explored in the literature. Layout is defined as the arrangement of stiffeners on the base shell. The Solid Isotropic Method with Penalization (SIMP) method (Bendsøe 1989) is applied for optimizing the layout of stiffeners or reinforcements in several topology optimization cases, for example, in maximizing the overall stiffness (Afonso et al. 2000) or eigenfrequencies (Du and Olhoff 2007). A different approach is the thickness optimization of the finite elements of the base shell. In this approach, the areas that have a higher 1 3 value than the set thickness threshold can be considered as potential stiffener locations (Lam and Santhikumar 2003;Khosravi et al. 2007). Also, the level-set method has been used to identify the optimal stiffener regions (Gomes and Suleman 2008). It should be noted that all aforementioned methods are suitable for identifying regions where stiffeners could be placed potentially. However, neither information about the sizes nor the topology of the stiffeners is retrieved. Usually, after interpretation of the results, a separate sizing optimization needs to be performed (Afonso et al. 2000).
For the simultaneous optimization of the stiffener layout and their size, three methods have been proposed in the literature. The first is the group of biologically inspired methods, such as the Adaptive Growth Method Yamazaki 2004, 2005;Dong et al. 2019Dong et al. , 2018Shen et al. 2019). To overcome the shortcoming of an empirical approach with user-defined parameters, which is unable to handle multiobjective cases, the Adaptive Growth Method was reformulated in terms of analytical rules that cover the morphogenesis of the growth of leaf veins in nature (Li et al. 2013). The second method is the method of Moving Morphable Components, applied to optimize stiffeners on a plate for maximum stiffness, or minimum compliance, subject to a volume and buckling constraint (Zhang et al. 2018). The last method is the Ground Structure Approach (GSA) (Dorn et al. 1964). This method was, for example, applied to the optimal panel placement in an airplane wingbox (Sleesongsom and Bureerat 2013; Yang et al. 2016).
The topology optimization problem for individual stiffeners, where the stiffener layout is pre-assumed, has been investigated in previous research, in particular, for applications to an airplane wingbox. Here, a ground structure of stiffeners is assumed and the topologies of the stiffeners are optimized using the SIMP method. The minimization of the compliance with a volume constraint was performed (Krog et al. 2004). Also, different constraints such as lift, drag, and stress for minimizing the mass were considered (Maute and Allen 2004). Optimization for a flutter and compliance objective under a weight constraint was performed in Stanford and Dunning (2015). Most recently, for this wingbox application, an optimization of the individual stiffener topology for minimization of the mass under buckling and stress constraints was reported (Stanford 2018). A more general application to stiffened panels was considered for a buckling objective with a volume constraint (Stanford et al. 2014). Recently, a level-set approach was published for stiffened plates with a pre-assumed stiffener layout. The topology of the individual stiffeners was optimized for a buckling objective under a mass constraint (Townsend and Kim 2019).
Modular structures presented in the previous research mainly focus on topological periodicity. In this setting, the design domain is divided into sub-domains which are constrained to be topologically identical. As such, a single module consisting of a ground structure of trusses is optimized for minimal weight, while subject to a fixed number of trusses (Beghini et al. 2014). For two-dimensional (2D) continua, a repeated module was incorporated by the use of a simple mapping technique, which separates the design variables of a module unit and the sub-domains in the global density field. The mapping from the design variables in the module unit is carried out one-to-one to the corresponding element material density values in the sub-domains (Almeida et al. 2010). There has been limited work available to date. An optimization of a stiffness criteria using the bidirectional evolutionary structural optimization (BESO) is performed (Huang and Xie 2008). This approach was later extended to natural frequencies (Zuo andXie 2011), conductivity (Chen et al. 2010), and a multi-objective formulation (Thomas et al. 2020). In the latter work, also a SIMP approach is utilized and connectivity between modules is taken into account. To improve computational efficiency for the design of a module, static condensation is utilized in combination with a level-set approach for a stiffness objective (Fu et al. 2019).
Although the aforementioned methods are able to design a structure consisting of a simple module repeated several times in the global domain, they suffer from a common limitation. Namely, the designs converge towards solutions with compromised structural performance (Huang and Xie 2008;Zhang and Sun 2006). The cause lies within the topological periodicity. The topology of the module is influenced most by the region with the highest compliance. The resulting module design is used at different locations in the structure, therefore not leading to the most optimal solution for these regions (Tugilimana et al. 2019).
This shortcoming can be addressed by two approaches: (i) by defining additional module properties as design variables or (ii) by allowing more modules within the structure. Both approaches extend the solution space. The first approach, extending the solution space, was considered by allowing for rotation of a module. Allowing for rotations resulted in improved structural performance because rotation of the modules modifies the material distribution in the structure locally (Tugilimana et al. 2017). Also, in a 2D continuum setting, the one-to-one mapping of a module unit to the global domain is extended by allowing the module unit to resize (Stromberg et al. 2011).
In the second approach, more than one module is allowed within the structure. The optimization problem is therefore redefined as the search for several module topologies and the distribution of these in the domain. This has been incorporated for truss structures based on the ground structure approach. Moreover, the modules were also allowed to rotate (Tugilimana et al. 2019). For a 2D continuum, the definition of a mapping between the design variables of a module unit and the corresponding sub-domains was extended to enable simultaneous optimization of multiple module unit topologies and their layout in the sub-domains (Higginson et al. 2020). The mapping is based on a weighted sum, allowing for the choice of one unique module unit in a sub-domain. The resulting topology optimization framework for modular structures can be combined with gradient-based optimization.
The aforementioned state of the art emphasizes the potential of structural optimization to enhance the conceptual design of structures. However, it is observed that research is mainly focusing on optimizing one of the following three aspects: (i) the stiffener layout, (ii) the individual stiffener topology, or (iii) truss-based or 2D continuum modules. Therefore, the goal of this paper is to develop an optimization method that simultaneously optimizes the modular stiffener components including their topology and layout on a base shell or plate.
The proposed method relies on the combination and extension of two existing methods: a ground structure approach which is combined with the topology optimization framework for 2D continuum modular structures (Higginson et al. 2020). The present work will focus on maximizing the overall stiffness of the structure, stated as minimizing the compliance, while subject to a prescribed volume. However, it is emphasized that the proposed method can easily be extended to other settings. The presented examples, from an engineering point of view, shall also be extended with at least vibration and buckling considerations. The main idea will be described on the basis of a simple example, consisting of a plate with stiffeners, see Fig. 1a. On the base plate, a ground structure of stiffeners is placed. For this example, the ground structure is generated by specifying a uniform grid on the base plate. In Fig. 1a, a ground structure for the stiffeners has been presented consisting of two stiffener types. The aim of the proposed method is to specify a fixed, but limited number of modules within these stiffener types and to find their optimal topology. The topologies of these modules can range from empty, called void, to fully present. These modules can be repeatedly used in the ground structure. The layout of the modules in the ground structure is simultaneously optimized with the topology. As such, a layout in the ground structure arises which only consists of a limited number of modules, as illustrated in Fig. 1b. The final topology optimization will be based on a SIMP formulation.
This paper is organized as follows: in Section 2, the detailed description of the proposed method is provided. In Section 3, the method is applied to several practical cases. Conclusions and recommendations are provided in Section 4.

Modularity in the ground structure
The proposed optimization method is based on a combination of a ground structure approach with the concept of material density topology optimization for modular structures. The main idea, as introduced in the Introduction, will be further specified on the basis of an example, consisting of a plate with stiffeners, see Fig. 2a. On a base plate, a ground (a) On top of a base plate a ground structure of stiffeners is generated. This ground structure consist out of two types of stiffeners which can be distinguished by their type letter.

(b)
The proposed method aims to find the optimal topologies of a fixed and limited number of module stiffeners for every type of stiffener. In this case, the number of modules for each stiffener type is limited to 3. These modules can be repeatedly used in the ground structure. The layout of the modules in the ground structure is optimized simultaneously with the topology. As such, a layout in the ground structure arises which only consists of the specified modules.

Fig. 1
Overview of the proposed method explained in two parts. In a an initial ground structure is presented, which is the basis for the optimization. A result of the layout of the modules within the ground structure and the topologies of the modules is shown in b structure of stiffener components with parents and children, each occupying a domain s , is generated. For clarity, in this case, the topology of the base plate will not be optimized and is therefore assigned as non-design domain n . A parent-children scheme is introduced in this section. At the generation of the stiffeners in the ground structure, one or multiple parent stiffener can be specified. In Fig. 2a, these parent stiffeners are the stiffener domains s=1 and s=6 . The parents can be identified by their type letter, in Fig. 2a, Type A and B, respectively. A mesh is generated for each parent stiffener, as highlighted in Fig. 2a. Since the topology optimization is SIMP based, a material density, is assigned per finite element e, e . The parents are copied one-to-one in the ground structure to form the so-called children for each type. In Fig. 2a, for Type A these children are stiffener domains s=25 and for Type B, stiffeners s=7−10 . The result is an initial mesh consisting of a base plate and a ground structure of parent stiffeners with their children.
If this mesh is subjected to a standard topology optimization, a unique topology is allowed to arise for every stiffener. As stated before, it is beneficial to limit the topologies of the parents and their according children to a fixed number of modules. Therefore, module templates are introduced. Module templates are an integer number T s , of one-to-one copies of the parent stiffener. As such, also these module templates have identical mesh topologies and inherent material densities as the according parents, see Fig. 2b. The material densities in the templates are defined using the finite element number in the template, d, and the module template number, t, denoted by t,d , see Fig. 2b. These template densities, t,d , are considered as the primary design variables and form the basis for the topology optimization. More details on the topology optimization are provided in Section 2.2.
The material densities of the mesh, e , will be determined by a mapping between the material densities of the templates t,d by the use of their according weight factors. The use of weight factors is inspired by the field of Discrete Material Optimization (DMO) (Higginson et al. 2020;Stegmann and Lund 2005). Here, a multi-material optimization is commonly described by an element constitutive matrix defined as a weighted sum of predefined potential materials. If a weight factor is a value around 1, the corresponding material is present in the element, if the value is around 0, the material is absent. This idea is used with the module templates and the parent-children scheme. A number of weight factors w s,t are assigned to each stiffener of a certain type. The number of weight factors is equal to the number of templates introduced for this type T s . For the example in Fig. 2b, 2 templates per type have been introduced. The weight factors, w s,t , denote the presence or absence of template t in stiffener domain s . Details on the mapping are provided in Section 2.3. By optimizing the weight factors simultaneously with the material densities of the templates, a layout (a) On the blue base plate, a ground structure of stiffener components with parents and children, each occupying a domain Ωs, is generated. For each stiffener type, here denoted by Type A and B, the generation of the stiffeners in the ground structure is based on a parent stiffener, in this example these are domains Ωs=1 and Ωs=6. As magnified, a mesh is generated for each parent stiffener respectively. Every finite element, e, is assigned a material density, ρe. The meshed parents are copied one-to-one to form the so-called children. For Type A these children are stiffener domains Ωs=2−5 and for Type B, stiffeners Ωs=7−10.
(b) To limit the topologies of the parents and their according children to a fixed number of modules, templates are introduced. Templates are an integer number of one-to-one copies of the parent stiffener. As such, also these templates have identical mesh topologies and material densities as the according parents. The material densities of the templates are defined in terms of the finite element number in the template, d, and the template number t, ρ t,d . These templates can be repeatedly used in the ground structure. The material densities of the mesh in the ground structure, ρe, will be determined by a mapping between the material densities of the templates ρ t,d depending on weight factors. The weight factors, ws,t, denote the presence or absence of template t in stiffener domain Ωs.

Fig. 2
In a the initialization of a ground structure is presented. In b templates are introduced in the ground structure arises, which only consists of the specified templates.

Topology optimization using the Solid Isotropic Method with Penalization
The topology optimization in this work is based on the SIMP approach (Bendsøe 1989;Zhou and Rozvany 1992). Each finite element is assigned a continuous pseudo material density , which is used as a variable to interpolate the Young's modulus E. With an initial Young's modulus E 0 for a linear isotropic material with Poisson ratio , the relation is Hereafter, the pseudo material density will be referred to as material density or density. The material density has a penalty factor p ≥ 1 . This implies that whenever the penalty factor is greater than 1, densities which are intermediate values of the range [0, 1] are penalized. The value of the penalization factor p is gradually increased from 1 to 5 during the optimization process. This so-called continuation approach drives the design gradually to a more distinct 0-1 design (Eschenauer and Olhoff 2001). In Section 2.3, the conditions for this continuation are discussed in more detail.
In the numerical examples presented in Section 3, the base plates and shells are assigned as non-design domains, n . Excluding these domains from the topology optimization can be enforced by setting the material densities for these finite elements to 1. It should, however, be emphasized that this is a choice, and also these base plates and shells could be part of the topology optimization. The volume of the resulting design in the total domain is now represented by In the current formulation, the results of the SIMP approach are not only dependent on the penalty factor p, but also on the size and orientation of the mesh (Eschenauer and Olhoff 2001; Maute and Allen 2004). However, this drawback can be removed by the use of filtering. This will be discussed in Section 2.6.

Mapping of the module templates to the ground structure
The modularity concept is introduced by introducing an element-based mapping scheme using the module templates.
The formulation does not only allow for optimization of the topologies of the module templates, but also for the optimal layout of these within the ground structure (Higginson et al. 2020). (1) The material densities of the mesh, e , will be determined by a mapping between the material densities of the templates t,d and by the use of their according weight factors w s,t . The material density of an element e ( ) of a certain type can be obtained from the material densities of the corresponding element in the templates t,d by Higginson et al. (2020): where e denotes the finite element number in the global mesh and s the according stiffener domain in the global mesh. A number of templates, T s , are introduced for this type, which are numbered by t. The finite elements in these templates are numbered by d. In this formulation, q > 1 denotes a penalty for the weight factors. If the factor is greater than one, intermediate values of the weight factor will be penalized. This scheme is similar to the penalty factor as used in SIMP, see (1). Therefore, the same continuation scheme is applied. In this work, the continuation is performed when three conditions are met. Firstly, the condition that the absolute change in the objective in two successive iterations is less than 0.1. Secondly, the designs should satisfy all the constraints during these two iterations. Finally, the last continuation should be more than 20 iterations ago.
An example of the mapping is provided for the problem in Fig. 3. The material density of dashed element e = 1 in stiffener domain s=1 is calculated. The stiffener is of Type A, so according to the mapping in (3), the material density of the element is determined by the dashed template elements, see Fig. 3: Fig. 3 The mapping of the topology optimized module templates to the mesh is illustrated through an example as earlier presented in Fig. 2. For the dashed element in the mesh, the weight factor w s=1,t=1 , corresponding to template t = 1 of Type A, is around 1. As such, the material density of this template is mapped to the mesh, resulting in a fully present element Some notes on the mapping should be made. First of all, it becomes clear from this formulation that in case a weight factor in the numerator is valued 0 or 1, the mapping does not provide a finite solution. Therefore, the weight factor should be limited to the range w s,t ∈(0, 1) . Secondly, as also noted by (Stegmann and Lund 2005), it could be observed that the term (1 − w s,t≠j ) q forces the design to a 0-1 solution for the templates, since an increase of one weight variable, automatically means a decrease in all other weights. Therefore, the converged values for the weight factors should denote a value around 1 if a template is present at a certain stiffener domain and 0 if it is not. Finally, the normalization term ensures that the overall mapping sums to unity for each stiffener domain.

Finite element analysis
The finite element analysis has multiple functions in the optimization. Primarily, it is used to model the physics. By implementing the boundary conditions, such as loadings and supports and successively evaluating the model, a response in terms of the global nodal degrees of freedom, denoted by , is retrieved. Using this information, the second function can be fulfilled: calculation of the objective and volume constraint. For this work, the objective is to maximize the overall stiffness of the structure. This can be stated as minimizing the compliance c, defined as where is the global stiffness matrix, which is a function of the element densities introduced by the SIMP approach. The well-known finite element equilibrium equation is where, denotes the nodal loads. In order to avoid singularities in the global stiffness matrix, the lower bound of the material density is chosen to be a small value min . In this work, this value is 5 ⋅ 10 −3 .
Lastly, the finite element analysis provides the sensitivities from the objective and the constraint with respect to (w.r.t.) the material density of the element. This will be further discussed in Section 2.8.
In this work, a triangular 6 node, 12 degrees of freedom, shell element will be used (Van Keulen and Booij 1996).

Problem definition
The search for the optimal structure with minimal compliance, while subject to (s.t.) static equilibrium, a maximum volume, a non-design domain n , and stiffener domains s in the global coordinate system using the SIMP formulation as in (1) and modular template mapping in (3), is stated as It becomes clear that the standard SIMP topology optimization problem, with design variables e ( ) is reformulated in terms of weight factors between the stiffener domains and templates w s,t and the element material densities of the templates t,d .

Density filtering per stiffener domain
Filtering is necessary in the SIMP approach to avoid solutions which are dependent on the mesh or provide checkerboards. One of these filters is the density filter, where the element densities are adjusted based on the values of the neighboring elements (Andreassen et al. 2011). An illustration is given in Fig. 4. A filter radius, dependent on a scalar value times the element size, is considered. For this work, this scalar value is set to 1. All the material densities of the elements, with their center-to-center distance Δ i within the relative radius r rel. , are weighted depending on this distance. Due to the triangular shape of the elements used in this work, the filter will include sufficient elements to overcome the checkerboard and mesh dependency issues, while allowing for small feature sizes in the topologies.
However, the standard density filter needs to be adjusted on three aspects to avoid mixing of the domains. First of all, to avoid material density values below 1 in the non-design areas n . Secondly, to prevent material density transfer from the non-design areas n to the stiffener domains s . Thirdly, to avoid transfer of material density through the edges of two stiffener domains s . Especially, if somewhere in the structure a full solid template is used, next to a neighboring void template, this would result in material density transfer between these and therefore destroy the clear boundary between solid and void. These three adjustments can be imposed by only taking elements in the same stiffener domain s , into account. The density filter is formulated for every element e in a mesh consisting of a total of N elements. The density filter takes into account the set of elements i, for which (i) the center-to-center distance i is smaller than the filter radius r rel. and (ii) is part of the same stiffener domain s . In equation form, this is stated as (Andreassen et al. 2011): where H i denotes the filter weight factor. This weight factor is defined as the maximum distance between the centers of the elements in s (Andreassen et al. 2011): In this work, the filter is applied at the design variables before entering the finite element analysis. The inverse of the filter is applied after the calculation of the sensitivities w.r.t. the objective and constraint functions, which will be discussed in the next section.

Gradient-based optimization
In this work, the optimization solver will be gradientbased using the Method of Moving Asymptotes (MMA) (Svanberg 1987). MMA is often used within topology optimization problems and has proved to be reliable in combination with multiple complex constraints (Koppen 2017).
In order for the optimizer to work properly, a scaling of the sensitivities might be required. The objective and constraint values should be scaled between 1 and 100 (Svanberg 1987). In this work, these values and their corresponding sensitivities are scaled by a constant to meet this criterion. The optimization procedure is terminated when the continuation of the penalty factors has reached 5 and a maximum number of iterations is reached or when the relative objective change is smaller than a prescribed amount. In this work, these are set to 400 iterations and 1 ⋅ 10 −7 .

Sensitivity analysis
The sensitivities of the design variables are calculated using the chain rule, since the stiffener domains get their designs from the according templates. Every material element density e is a function of the template material element density t,d and the according weight factor w s,t through the mapping as described in (3). Therefore, the sensitivity of the objective w.r.t. the element material density of a template can be written as where S t denotes the subset of elements e, which retrieve their material densities from the template. For example, if template t = 1 of Type A is used n times in the domain, then S t contains n values. The sensitivity of the objective w.r.t. the element material density c e ( ) is calculated using the adjoint formulation (Bendsøe and Sigmund 2003).
Due to the mapping as described in (3), the sensitivity of each element material density w.r.t. the template element material density of the same type is The sensitivity of the objective w.r.t. the weight factor is determined by summing each contribution of the corresponding template element material density d: The sensitivity of each element material density w.r.t. the weight factor can also be determined by taking the derivative of the mapping as described in (3). This is done in a similar fashion as in (11) and is therefore omitted.

Fig. 4
A detailed view on the filtering at the boundary of three domains. Here, stiffener domains s=1 and s=6 and the non-design domain n meet each other. The filter takes into account all the elements that are, with their center-to-center element distance Δ i , in the filter region with relative radius r rel. . The relative filter radius is dependent on a scalar times the element size. In order to prevent mixing of the material element densities over the boundaries of the domains, the filter is restricted to only take into account elements that are within the same domain. This is indicated by the green circle segment of the filter. The red circle segment indicates elements which are not part of the same domain The sensitivities of the volume constraint w.r.t. the template element material densities of the same type and the weight factors can be determined similarly: where the sensitivity of the volume constraint w.r.t. the element material density V e ( ) can be calculated by taking the derivative of the formulation described in (2).

Initial conditions
The initial conditions for both the weight factors and the template material densities are given for the first iteration. As stated in Section 2.3, the weight factor should be in the range (0, 1) to provide a finite solution. Therefore, the initial values for the weight factors are selected as where T s is the amount of introduced templates of this type.
In topology optimization, it is quite common to take 1 or the volume fraction as initial value for the element density variables. In the latter, the volume constraint is directly satisfied. However, it should be noted that if all element material densities have the same value, the sensitivities w.r.t. the weight factors will be equal for all the templates. As an example, consider Fig. 2. The initial value for the weight factors, according to (15), is 0.5 for both types. If the initial material density of the elements is also 0.5, the sensitivities of the objective and constraint w.r.t. the weight factors, respectively (12) and (14), would lead to sensitivity values equal to 0. As such, due to the gradient-based optimizer, no further changes in the weight variables or template element material densities would be posed. The result would be a solution, where every stiffener will have an equal contribution of the two templates with 0.5 as element material density, which removes the introduced modularity.
To give the sensitivity values a different direction and magnitude, the initial values of the element densities are slightly varied. This enables the optimizer to converge to distinct topologies and layout in the ground structure.
here t is the template number and T s the number of introduced templates of this type.

Numerical examples
This section discusses three numerical examples. The first two are examples of stiffened plates. Here, different initial ground structures, load cases, and parent-child schemes are used. The third example is a stiffened shell representing a fuselage of an airplane.

Simply stiffened plate
In Fig. 5, a simply stiffened plate is shown. The edges of the stiffeners and base plate are fully clamped. At the center of the base shell, there is a concentrated force F 1 . The values for the parameters are shown in Table 1. The parent-child scheme and the according stiffener domains are already introduced in Section 2.3 and shown in Fig. 2. The optimization problem formulation in (7) is utilized with an upper value of half the initial volume of the stiffeners for the volume constraint, V max = 0.5 ⋅ V init.stiff. . The influence of the number of templates per type is investigated. The example is performed for three cases, where the number of templates per stiffener type is varied from 1 to 3. Since the problem is symmetric along the x 1 − and x 2 -axis, it is hypothesized that the resulting stiffener layout and template topologies should be symmetric. The presented mesh has 184471 nodal Fig. 5 Geometry, loading, and boundary conditions for simply stiffened plate. The values of the parameters are given in Table 1 points and 91800 elements. The meshes used are added as Supplementary data.

Optimization results and discussion
The resulting template topologies of the three cases are shown in Fig. 6 and their layouts in the stiffener domains are shown in Table 2. The results and convergence data are added as Supplementary data. The observation is made, that the resulting template topologies and their layout in the stiffener domains are always symmetric. The case with 1 template has the highest compliance value and therefore the worst mechanical performance. As already stated in the Introduction, due to the topological periodicity, the template is influenced the most by the highest loaded region. In this case, these regions are the stiffener domains s=3,8 , which carry a majority of the concentrated load. The volume constraint does not allow for a complete solid topology of the templates and therefore an inferior solution arises. In the case of 2 and 3 templates per type, one complete solid template arises, which carries the major part of the concentrated load. The compliance values of these cases are therefore very comparable, with a minor improvement of the objective function for the 3 templates per type case. The compliance value of the latter two cases is around 200% lower than for the 1 template per type case. Since the plate consists of 5 stiffeners per type, it could be expected that up to 5 unique templates could be defined with even lower compliance than for the case of 3 templates. However, for the cases of 4 and 5 templates, similar topologies and layouts as for the 3 templates case are found due to the symmetry. This is characterized in the results by duplicate or not used templates, therefore the results of these cases are not presented. The results and convergence data are added as Supplementary data.

Orthogonally stiffened plate
A base plate with an orthogonally ground structure of stiffeners is considered. The stiffener domains are all based on one parent, as shown in Fig. 7a, and therefore only one template type is defined. The geometric features including the distributed load, concentrated force, and boundary conditions are shown in Fig. 7b. The entire base plate with all the stiffeners is modeled. However, it should be noted that the problem is symmetric. Therefore, it is hypothesized that the resulting topology should also be symmetric along the x 1 and x 2 planes. The optimization problem formulation as stated in (7) is used. The upper value of the volume  Fig. 6 An overview of the topologies of the templates for the simply stiffened plate example. In a the layout of the templates in the stiffener domains is shown for the case of 3 templates per type. The template topologies for this case are shown in b. For the case of 2 templates and 1 template per type, the resulting topologies are shown in c and d, respectively. Their layouts in the domains are given by Table 2 1 3 constraint is set to one-third of the initial volume of the stiffeners, V max = 0.33 ⋅ V init.stiff. . The number of templates per type is varied, with the hypothesis that an increased number of templates results in a lower compliance value, since the design space is increased. The presented mesh has 237871 nodal point and 118800 elements. The meshes are added as Supplementary data.

Optimization results and discussion
The resulting template topologies for the cases of 1 to 7 templates per type are shown in Fig. 8. For all cases, the resulting layout is symmetric along the x 1 -and x 2 -axis, therefore to keep it clear, only the results of one-quarter of the plate are shown. In this quarter plate, symmetry is also observed. The template layout is given for two opposite stiffener domains in Table 3. The lowest, and very comparable, compliance values are retrieved for the case of 5, 6, and 7 templates per type. For the case of 7 templates per type, the resulting layout is shown in Fig. 8a. It is noted that a major part of the volume is assigned to the stiffener domains s=7−9 and s=16−18 . These domains represent the stiffeners crossing the center of the base plate and carrying the major contribution of the concentrated force F 2 . Since in the center of the base shell the sum of the deflection due to the concentrated load and pressure load is expected to be the largest, it is reasonable that most of the volume should be used for these stiffener domains. This observation can be extended to all other cases, where the templates with the most volume are used at these stiffener domains as well.
Increasing the number of templates resulted in a lower compliance value, as can be observed in Fig. 8. Since the quarter plate consists of 9 unique stiffener domain pairs, it could be expected that up to 9 unique templates could be defined with even lower compliance than for the case of 7 templates. However, for the cases of 8 and 9 templates, similar topologies and layouts as for the 7 templates case are found. This is characterized in the results by duplicate or not used templates, therefore, the results of these cases are not presented. The results and convergence data are added as Supplementary data.

Orthogonally stiffened shell: airplane fuselage
A practical example inspired on the top middle part of the fuselage of an airplane is considered, especially, the critical loadings during a 2.5G pull-up maneuver with a pressurized cabin at cruise height (Şen 2010). The geometry, boundary conditions, and loadings are shown in Fig. 9. The parameters used are denoted in Table 4. It should be noted that for correct boundary conditions, the slider at the lower side of the base shell along the x 1 -axis should be in radial direction as is the case for the upper side. In this case this is not (a) A quarter of the orthogonally stiffened plate is shown, since the result is symmetric along the x 1 and x 2 axis. The parent is the stiffener domain Ω s=1 , shortly denoted as Ω 1 . All the other stiffener domains are children from this parent. The base plate is assigned as non-design domain Ω n .
(b) A orthogonally stiffened plate, subjected to a distributed load p, concentrated load F 2 and fully clamped on the edges of the base shell and stiffeners domains. The values of the parameters are given in Table 1. implemented, since the goal of this example is to show the application of the developed method to shell structures. The optimization problem formulation as stated in (7) is used, with an upper value of the volume constraint set to two-third of the initial volume of the stiffeners, V max = 0.667 ⋅ V init.stiff. . The only domains that are subject to the optimization are the webs of the stiffeners. The top of a stiffener and the base shell are assigned as non-design domain. The example is performed for three cases, where the number of templates per stiffener type is varied from 1 to 3. The presented mesh has 212193 nodes and 105910 elements. The meshes are added as Supplementary data.

Optimization results and discussion
The results of three cases are shown, where 1 up to 3 templates per type are used. For all cases, results and convergence data are added as Supplementary data. For the case of 1 template per type, the topologies are shown in Fig. 10d and the layout in the stiffener domains is given in Table 5. For the stiffeners of Type A, the first observation can be done that a fully present stiffener is arising. This is probably due to the fact that the p 2.5G load is applied at the side of the structure on the stiffeners of Type A at domains s=16→20 .
Removing these stiffeners will result in a structure, which is not able to transfer the load and is therefore avoided by the optimizer. A second observation can be performed from the stiffeners of Type B. Due to the volume constraint, material must be removed. A major part of this material is removed just below the stiffener top. This makes sense, since carrying the p press. by the structure can be done most effectively by placing material as far from the origin, and therefore as close to the base shell as possible.
In the case of 2 templates per type, the template topologies are shown in Fig. 10c and the layout in the stiffener domains is given in Table 5. The observations of the previous case can be extended. A void template t=1 for Type A arises, which is used mainly at domains s=6→9 and s=11→14 . Since the axial load p 2.5G is dominant in magnitude, these tangential stiffeners are removed. A similar observation can be done for the almost void stiffener template t=1 of Type B that arises. This stiffener template is used at the boundaries of the structure, where the boundary conditions provide stiffness. Moreover, as observed previously, more material is placed as close to the base shell as possible in stiffener template t=1 of Type B.
For the third case of 3 templates per type, the template topologies are shown in Fig. 10b and the layout in the stiffener domains is given in Table 5 and shown in Fig. 10a. It can be observed for both stiffener types that an additional template arises, which places material as close to the base shell as possible. For the stiffener template t=2 of Type B it is observed that it is used at the stiffener domains s=2→7,3→8,4→9 . Here, the boundary condition is applied at the base shell at domains s=1→5 . The majority of the p 2.5G load is carried by the base shell, and therefore, some more material is removed below the stiffener top as compared to template t=3 . Overall, it can be observed that the compliance decreases with the introduction of more templates.

Conclusions and recommendations
Stiffened shells are widely used in engineering structures, but their performance is highly influenced by the topology of the stiffeners and their layout on the base shell. Moreover, it is beneficial to design structures with modules, since it allows for increased and cheaper quality control, more accessible mass production, and therewith reduction of costs. The method proposed in this work allows for the simultaneous optimization of the topology of the modular stiffeners and their layout on a base shell or plate.
It can be briefly concluded that the proposed method -enables simultaneous optimization of the stiffener topology and layout of stiffeners on shells and plates; -incorporates a fixed, but limited number of integer modules in a three-dimensional structure; -reduces the number of design variables; -prevents mixing of the boundaries of the domains through the adjusted density filter; (a) Stiffener layout for the case with 7 templates per type.  Fig. 8 An overview of the topologies of the templates is shown for the orthogonally stiffened plate. In a the layout of the templates in the stiffener domains is shown for the case of 7 templates per type. The template topologies for this case are shown in b. For the case reaching from 1 to 6 templates per type, the resulting topologies are shown in c-h. Their layout in the domains is given in Table 3 Fig. 9 Geometry, loading, and boundary conditions for the orthogonally stiffened shell inspired on an airplane fuselage. The values of the parameters are given in Table 4. The stiffeners are located at the inner side of the base shell segment, at the dashed lines. The parent for Type A is the stiffener domain in tangential direction between Nodes 1 and 2, denoted by s=1→2 . For Type B, the parent is the domain in axial direction between Nodes 1 and 6, denoted by s=1→6 -gradually drives the module templates to topologies with distinct solid/void boundaries and a clear layout in the ground structure because of the continuation scheme; and -allows conceptual designs for stiffened plates and shells with different number and formats of the module templates.
The proposed method is applied to several stiffened plate examples and a practical stiffened shell problem. For all examples, the method has shown to converge to designs with distinct solid/void boundaries in the module topologies and a clear layout in the ground structure. These designs can be reasoned logically given the boundary conditions and loads. For example, the symmetric stiffened plate examples resulted in symmetric designs. Moreover, it was shown that allowing for more templates increased the design space and therefore resulted in lower values for the compliance. Even for a limited number of templates and therefore a reduced number of design variables, the method is also able to generate innovative designs. Therefore, the proposed method is a design tool that can be utilized in the conceptual design phase of structures with stiffeners. The recommendations are separated into two parts: the current implementation and future applications or extensions of the method. For the first part, it is recommended to develop a more thorough understanding of the influence of the continuation scheme on the final results. The conditions for increasing the penalty factors and their influence on the final design should be more thoroughly investigated. The same holds for the filtering. In this work, an adjusted density filter with a constant relative filter radius is adopted. The filter radius implicitly describes the minimal feature size that can arise in the topology. Therefore, the influence of the filter radius has to be taken into account in more detail. However, the density filtering also opens up new possibilities. For example, the filter radius could be adjusted per module template to provide a minimal feature size per module.
For the second part, future applications and extensions of the method, the importance of different objective functions is recognized. As already stated in the Introduction, thinwalled structures are sensitive to vibrations. Therefore, it is recommended to include the dynamics in the objective functions. As also stated in the Introduction, the structural performance of modular structures can be enhanced by allowing additional module properties as design variables or by introducing more modules to the structure. This work only focuses on the latter and therefore the inclusion of additional module properties such as the rotation can be incorporated in future work. In this work, the method is mainly applied to stiffener domains. It should be noted that the method is defined very generally. This opens up possibilities to apply the method to domains of different shapes to further extend the range of applications. Moreover, the generality of this method allows for different topology description methods. In this work, the material density-based topology optimization using SIMP was used, although the method can also be combined with other topology description methods, such as level-set.  Fig. 10 An overview of the topologies of the templates for the orthogonally stiffened shell. In a the layout of the templates in the stiffener domains is shown for the case of 3 templates per type. The template topologies for this case are shown in b. For the case reaching from 2 and 1 templates per type, the resulting topologies are shown in c, d, respectively. Their layout in the domains is given in Table 5 Supplementary Information The online version of this article contains supplementary material available https:// doi. org/ 10. 1007/ s00158-021-03081-0.