Multi-Domain and Multi-Material Topology Optimization in Design and Strengthening of Innovative Sustainable Structures

Abstract

Expectations and challenges of modern sustainable engineering and architecture stimulate intensive development of structural analysis and design techniques. Designing durable, light and eco-friendly constructions starts at the conceptual stage, where new efficient design and optimization tools need to be implemented. Innovative methods, like topology optimization, become more often a daily practice of engineers and architects in the process of solving more and more demanding up-to-date engineering problems efficiently. Topology optimization is a dynamically developing research area with numerous applications to many research and engineering fields, ranging from the mechanical industry, through civil engineering to architecture. The motivation behind the present study is to make an attempt to broaden the area of topology optimization applications by presenting an original approach regarding the implementation of the multi-domain and multi-material topology optimization to the design and the strengthening/retrofitting of structures. Moreover, the implementation of the design-dependent self-weight loading into the design model is taken into account as a significantly important issue, since it influences the final results of the topology optimization process, especially when considering massive engineering structures. As an optimization tool, the original efficient heuristic algorithm based on Cellular Automata concept is utilized.

1. Introduction

Observed for decades, the continuous and intensive development of industry and construction has led to both positive and negative changes in the many areas of human life. Nowadays, the efforts of scientists, engineers, architects and planners are focused on building innovative sustainable architecture restoring the balance between the products of human hands and nature as nature is thought to have created the best environment for humans. To achieve the goal, architecture, urban planning, civil engineering, material science, computer science and mechanical engineering cooperate to make good designs that can improve the quality of life. As a result of this cooperation, new innovative and efficient design methods are developed, which can combine aesthetics and ecology with reasonable engineering requirements and limitations of economy.

Among innovative design tools, topology optimization took its place as an important approach, soon after it was introduced in [1]. This powerful and versatile approach allows original designs to be created by an optimal distribution of material within the design space. At the same time, over the last decades, one could observe the implementation of various specific topology optimization methods, ranging from gradient-based approaches, e.g., [1,2,3,4], where mathematical models are derived to calculate the sensitivities of design variables, to non-gradient-based ones, where material is redistributed using various, in many cases heuristic, techniques. In what follows, generation of optimal topologies involves, among others, Evolutionary Structural Optimisation (ESO) and Bi-Directional Evolutionary Structural Optimization (BESO), e.g., [5,6]; genetic algorithms, e.g., [7,8]; other biologically inspired algorithms, e.g., [9,10]; the material cloud method, e.g., [11]; spline-based topology optimization, e.g., [12]; the level set method, e.g., [13,14,15,16]; proportional topology optimization [17]; or moving morphable components approach, e.g., [18,19].

The most recent development is observed when coupling the Solid Isotropic Material with Penalisation (SIMP) approach with the Non-Uniform Rational Basis Spline (NURBS) hyper-surfaces framework. The robustness and effectiveness of the SIMP method combined with the NURBS entities provide a pure geometrical (CAD-compatible) descriptor of the topology of the continuum. The isogeometric topology optimization method [20,21,22] has recently been used to solve problems involving requirements on structural displacements and eigenfrequencies, as well as in multi-scale topology optimization problems. A variant of this method has also been used for the multi-material topology optimization of continuum structures [23].

Many survey papers and books provide broad discussions on topology optimization concepts, e.g., [24,25,26] or [27]. At the same time, hundreds of papers present numerous results obtained for both simple models and complicated engineering structures. Nevertheless, although remarkable achievements have been already made, there is still room for further investigations.

Topology optimization has been broadly accepted by the design community, and literature presents various civil engineering and architectural applications like designing bridges, e.g., [28,29]; designing high-rise buildings, e.g., [30]; and others e.g., [31]. The discussion regards various issues including, for example, different types of constraints [32] and the optimization of bi-material structures [33]. The strengthening of civil structures suffering from the effects of material degradation, an extension of the standard area of application of topology optimization, was proposed as well. This very important and up-to-date topic of retrofitting heritage structures has been reported in the literature and [34,35] may serve here as examples. Topology optimization allows original designs to be obtained for civil engineering and architectural structures and can complement, for example, tensegrity designs [36,37].

This paper proposes an original insight into the development of design and strengthening of civil structures within innovative sustainable architecture by introducing a heuristic approach to multi-material topology optimization, which can be performed for specially arranged structural sub-domains for which different materials are implemented. Once the design domain has been divided into sub-regions, optimal topologies are generated under the restriction that the redistribution of material can be performed only within sub-domains selected for the employed materials. This approach is an alternative to the standard one, which is the redistribution of one or a few materials within the whole design domain.

The presented concept is completed by including self-weight loading, which makes obtained solutions more practical and realistic. Topology optimization of structures under design-dependent load is rather rarely considered in the literature. The first consideration of self-weight loading in optimal design can be found in [38], one of the first applications of topology optimization to structures under self-weight loading is discussed in [39], and some architectural applications can be also found in, e.g., [40]. The implementation of the self-weight loading has a significant influence on the final results of topology optimization process, which is very important, especially when dealing with the optimization of massive engineering structures, like, for example, bridges or carrying systems of tall buildings. The nature of the design-dependent loading is that the location, direction and magnitude of the load depend on the actual shape of the structure. Among examples, one can point out surface pressure caused by wind loads, e.g., [41]. Topology optimization of structures under design-dependent loads can be a challenging task because it cannot be treated as an extension of standard formulation of topology optimization under fixed loads and therefore requires a modification of well-recognized approaches in order to eliminate issues like a non-monotonous behavior of the compliance or parasitic effect [1].

The outline of the paper is as follows. In Section 2, the topology optimization problem is formulated, and the concept of sub-domain-based optimization is described. The Cellular Automata topology generator is introduced in Section 3, where local update rules and the flexible approach are discussed. In Section 4, multi-domain and multi-material topology optimization of support and multi-layer structure under an external load and self-weight is performed, and obtained results are presented. The discussed approach is adapted also to strengthening and retrofitting selected structures in Section 5, where the numerical calculations are performed both for single- and multiple-load cases. The paper ends with concluding remarks.

2. Sub-Domain Based Multi-Material Topology Optimization

Within a design domain, the problem of structural topology optimization is to find the distribution of material giving the minimal value of the compliance of the structure c (Equation (1)). The design domain is discretized by finite elements. The available material volume fraction κ is defined and treated as a constraint in the optimization process (Equation (2)).

$$\mathrm{minimize} c\left(d\right)=U^{\mathrm{T}}KU={{\displaystyle \sum }}_{i=1}^N{d}_i^p{u}_i^T{k}_i{u}_i$$

(1)

$$\mathrm{subject} \mathrm{to} V=\kappa V_0$$

(2)

$$KU=F$$

(3)

$$0<d_{min}\le d_i \le 1.$$

(4)

In Equations (1)–(4), U and F are the global displacement and force vectors, K is the global stiffness matrix, u_i and k_i are the element displacement vector and stiffness matrix, respectively, and N is the number of elements, while d_i is a design variable that represents a relative density of the material assigned to each element. The constraints are imposed on the values of design variables, and d_min is a non-zero minimum relative density introduced to avoid singularity.

As the material representation, the SIMP (Solid Isotropic Material with Penalization) power law approach [42] is adapted. The elastic modulus E_i of each finite element is modeled as a function of the design variable d_i. In Equation (5), p (typically p = 3) penalizes intermediate densities and drives the design of black-and-white structures. E₀ and ρ₀ are the elastic modulus and density of a solid material, respectively.

$$E_i=d_i^p{E}_0, {\rho }_i=d_i{\rho }_0.$$

(5)

When dealing with the self-weight loading, the modification and extension of the standard SIMP method is necessary (see [43]), such that the solid material density ρ₀ of each element is defined according to Equation (6):

$${\rho }_i=d_i^p{\rho }_0.$$

(6)

It is proposed to extend the conventional concept of topology optimization by introducing the idea of multi-domain and multi-material topology optimization. A typical approach regards redistribution of one or several materials within a whole design domain, but this paper proposes a different approach. The design domain is divided into sub-regions for which different types of material are defined, and through this procedure, the multi-material structure is created. Following the preliminarily findings of [44], the aim of the present research is to find optimal topologies, under the restriction that redistribution of material can be performed only within sub-domains selected for employed materials. What is very important, in terms of practical applications, is that it is possible to impose separate constraints on volume fractions of each defined material. This approach produces different results when compared to the traditional problem, which is the redistribution of one or few materials within the whole design domain. Allowing for implementation of many materials during the topology optimization process may open new possibilities for improving existing solutions. In addition, including the self-weight loading makes the considered design problems more practical and realistic, although more complicated in terms of computation. The nature of this type of loading, namely, design-dependent loading, is that the location, direction and magnitude of the loads depend on the current shape of the structure in the current optimization step. Topology optimization of structures under design-dependent loads is rather rarely considered in the literature. This approach requires modifying known and widely recognized approaches due to the potential occurrence of such issues as a non-monotonous behavior of the compliance, parasitic effect or unconstrained character of the optimum, which means that the volume constraint can be inactive for optimal topology.

3. Cellular Automata Rules for Topology Optimization

The effectiveness of the topology optimization process is determined by selecting a proper method of numerical optimization. Heuristic optimization techniques are gaining popularity among researchers—see, e.g., [45,46]—because they are easy for numerical implementation, do not require gradient information, and one can easily combine this type of algorithm with any finite element structural analysis code. This practical aspect of engineering implementation of topology optimization techniques seems to be an important issue for contemporary design.

This paper utilizes an efficient heuristic approach based on the concept of Cellular Automata (CA). The main idea of CA is to replace a complex problem with a sequence of simple decision-making based on local information and local interaction between cells. The implementation of Cellular Automata requires decomposition of the design domain into a usually, but not always, e.g., [47], uniform lattice of cells that are equivalent to finite elements while performing analysis and topology optimization. It is assumed that the interaction between cells takes place only within the neighboring cells and is governed by local rules. The rules are identical for all cells and are applied simultaneously to each of them. Updating the state of cells can be performed based on already updated values found for cell neighbors (the Gauss-Seidel update scheme) or based on the states of the surrounding cells determined in the previous iteration (the Jacobi update scheme).

In this paper, a heuristic local update rule is proposed utilizing the Jacobi update scheme (see Equation (7)), where t is the current iteration number:

$$d_i^{\left(t+1\right)}=d_i^{\left(t\right)}+\Delta d_i^{\left(t\right)},$$

(7)

$$\Delta d_i^{\left(t\right)}=\frac{1}{M+1}\left[{{\displaystyle \sum }}_{k=1}^M{\alpha }_{ik}+{\alpha }_i\right]m.$$

(8)

In Equation (8), m denotes move limit (e.g., m = 0.2). The coefficients α_i for central cell i and α_ik for M neighboring cells are calculated based on the following rules [48]:

$${\alpha }_i=\left\{\begin{array}{c}-C_{\alpha } \mathrm{if} c_i \le c^{*}\\ \\ C_{\alpha } \mathrm{if} c_i>c^{*}\end{array}\right.; {\alpha }_{ik}=\left\{\begin{array}{c}-C_{\alpha } \mathrm{if} c_{ik} \le c^{*}\\ \\ C_{\alpha } \mathrm{if} c_{ik}>c^{*}\end{array}\right.,$$

(9)

where c_i (Equation (9)) is the value of the compliance of a central cell, and c_ik stands for the compliance of neighboring ones. The quantity C_α is a user-specified parameter [48].

In this paper, a modification of these rules is proposed. To make local rules more flexible, instead of one threshold value c*, two values c₁* and c₂* are introduced.

$${\alpha }_i=\left\{\begin{array}{c}-C_{\alpha } \mathrm{if} c_i \le c_1^{*}\\ \\ f\left(c_i\right) \mathrm{if} c_1^{*} \le c_i\le c_2^{*}\\ \\ C_{\alpha } \mathrm{if} c_i>c_2^{*}\end{array}\right.; {\alpha }_{ik}=\left\{\begin{array}{c}-C_{\alpha } \mathrm{if} c_{ik} \le c_1^{*}\\ \\ f\left(c_{ik}\right) \mathrm{if} c_1^{*} \le c_{ik}\le c_2^{*}\\ \\ C_{\alpha } \mathrm{if} c_{ik}>c_2^{*}\end{array}\right..$$

(10)

At the same time, the linear function in Equation (10) is built to fulfill the conditions f(c₁*) = −C_α and f(c₂*) = C_α; thus

$$f\left(c\right)=2C_{\alpha }\frac{c}{c_2^{*}-c_1^{*}}-C_{\alpha }\frac{c_2^{*}+c_1^{*}}{c_2^{*}-c_1^{*}}.$$

(11)

The values of the function f(c) (Equation (11)) refer to all compliances from the interval [c₁*, c₂*]. The width of the interval can be modified during the optimization process. This gives the function f(c) more flexibility, allowing for the choice between linear functions of different slopes, and the step (−C_α, C_α) one if c₂* − c₁* tends to zero. One can, for example, start with a wide interval, decreasing it as the generation of topologies proceeds, or the opposite: start with one threshold value c₁* = c₂* = c* and then increase the interval around c*. Implementing the adaptive technique by decreasing the interval [c₁*, c₂*] can, for example, speed up the process of the so-called grey elements elimination at the final stage of topology generation.

The numerical algorithm was built in order to implement the design rule proposed above. As for the optimization procedure, the sequential approach was adapted, meaning that for each iteration, the structural analysis performed for the optimized element was followed by a local updating process. Simultaneously, a global volume constraint can be applied for a specified volume fraction. The volume constraint is implemented in each iteration when local update rules are applied to all elements. In practice, the design variables multiplier is introduced, and then its value is sought to fulfill the volume constraint. As a result, the generated topologies preserve a specified volume fraction of a solid material during the optimization process. In this paper, the Finite Element Analysis tool (ANSYS) stands for the structural analysis module.

In the present study, the assumed change of the objective function value for subsequent iterations was implemented as the stopping criterion, but in general, the criterion can also be defined as performing a selected number of iterations.

It is worth noticing that the total volume constraint can be inactive while considering topology optimization including self-weight loading. The proper choice of c₁*, c₂* and C_α can eliminate this problem.

4. Multi-Domain and Multi-Material Topology Optimization

4.1. Example 1—Support Structure

4.1.1. Topology Optimization of Uni-Material Structure under External Load

As the first example, a plane elastic support structure bolted on two brackets presented in Figure 1 was chosen, assuming the symmetry of the structure. The structure is discretized with a regular lattice of tetragonal cells. The number of cells equals 8005. Firstly, the topology optimization using one material M₁ with the Young modulus equal to E₁ = 2 × 10¹¹ Pa, the Poisson ratio ν₁ = 0.3 and density ρ₁ = 8000 kg/m³ was performed. The volume fraction equal to 0.5 was implemented. The volume fraction was constant during the optimization process, and it is checked after each iteration (applies to all the examples considered in the paper). Due to symmetry, one half of the structure was considered.

For the first calculation, a concentrated force equal to P₁ = 1000 N was applied. The final topology with the resultant compliance equal to 4.866 × 10⁻⁶ Nm is presented in Figure 2, while the compliance of the initial structure with uniformly distributed material is equal to 8.482 × 10⁻⁶ Nm. At the same time, the maximal displacement of 0.243 × 10⁻⁷ m is reached for the final distribution of the available material, while for the initial structure the maximal displacement equals 0.424 × 10⁻⁷ m. It is important to note that the final topology for this case of loading does not depend on the value of applied loading.

To illustrate the effectiveness of the proposed approach, the computing time for the examples discussed in the paper was examined. As a metric of the algorithm runtime, a clock-time was used. Each iteration consisted of two stages: FEM analysis and design variable update. The time of the update step had a minor influence on the resultant time of the calculation process (it took less than 1 or 2 s, even for a large number of elements), while the structural analysis time depended highly on the mesh resolution (the number of elements used for the discretization of the structure). For all 2D examples considered in the paper, the 4-node finite element was used, and for 3D examples, the 8-node element. ANSYS Mechanical APDL was used as the analysis tool with AMD Phenom(tm) II X4 955 3.20 GHz processor and 8 GB of RAM (64-bit architecture). For the example presented in this section (8005 elements), 30 iterations were needed to reach a convergence, which took about 3 min (one half of the structure was investigated due to the symmetry). The iteration history showing the convergence of the objective function is provided in Figure 3.

In order to broaden the discussion regarding this example, the problem has been reformulated by replacing the volume constraint by the mass constraint. In what follows, the mass constraint µ = 0.5 has been implemented for the same geometry, material data and boundary conditions. As was expected for uni-material problem, the final values of compliance and maximal displacement are the same, since the volume and mass constraints are in this case equivalent.

The final topology, with resultant compliance and maximal displacement equal to 4.866 × 10⁻⁶ Nm and 0.243 × 10⁻⁷ m, respectively, is presented in Figure 4, whereas Figure 5 shows iterations history. It can be noted that less than 30 iterations were needed to reach a convergence.

4.1.2. Topology Optimization of Multi-Material Structure under External Load

For the next example, the two materials M₁ and M₂ with different material data were applied as shown in Figure 6.

Both materials are linear elastic and isotropic. Material 1 (M₁—drawn in Figure 6 in dark grey color) is characterized by a Young modulus equal to E₁=2 × 10¹¹ Pa, the Poisson ratio ν₁ = 0.3 and density ρ₁ = 8000 kg/m³; the second one—material 2 (M₂—drawn in Figure 6 in light grey color) is characterized by a Young modulus equal to E₂ = 2 × 10¹⁰ Pa, the Poisson ratio ν₂ = 0.3 and density ρ₂ = 2000 kg/m³. The ratios of the Young modulus and density values calculated for material M₁ and material M₂ are approximately the same as the ratios of these quantities for steel and concrete (compression). The same calculations as for the uni-material structure were performed, i.e., for a structure loaded only with a concentrated force equal to 1000 N. The final topology gained for this case is presented in Figure 7. The presented results show that multi-material topology optimization has to be treated as a separate/different approach, since it can give a new layout of material distribution for the optimized structure.

The final compliance of the optimized structure equals 0.310 × 10⁻² Nm, while the compliance of the initial structure with uniformly distributed material is equal to 0.539 × 10⁻² Nm. At the same time, the maximal displacement of 0.176 × 10⁻⁵ m is reached for the final distribution of the available material, while the maximal displacement of the initial structure is equal to 0.269 × 10⁻⁵ m.

For the example presented in Figure 7 (8005 elements), 33 iterations were needed to reach a convergence, which takes about 3 min (half of the structure was investigated due to the symmetry). The iteration history showing the convergence to the final solution is provided in Figure 8.

The presented example was remodeled into a three-dimensional structure by adding a third dimension—the thickness equal to 0.25 m. A quarter of the structure was investigated due to symmetry. All materials data as well as the boundary conditions are the same as for the 2D example presented in Figure 6. The same calculations were performed to illustrate the effectiveness of the algorithm. The value of the final compliance for the optimized structure equals 0.155 Nm (0.303 Nm for the initial structure), and the final maximal displacement equals 0.386 × 10⁻⁴ m (0.758 × 10⁻⁴ m for the initial structure). The final topology is presented in Figure 9.

For the example presented above (20,808 elements), 37 iterations were needed to reach a convergence, which takes about 12 min (a quarter of the structure was investigated due to the symmetry). The iteration history showing the convergence to the final solution is provided in Figure 10.

As for the previous example, the considerations regarding the multi-material 2D structure have been extended by replacing the volume constraint with the mass constraint. In what follows, the mass constraint µ = 0.5 was implemented for the same geometry, material data and boundary conditions. This time, since the multi-material structure is investigated, the constraint replacement influences the final results. The compliance of the optimized structure equals 0.276 × 10⁻² Nm, while the maximal displacement 0.138 × 10⁻⁵ m. The final topology is presented in Figure 11. The iteration history showing the convergence to the final solution is provided in Figure 12 (30 iterations were needed to get the solution). It is worth underlining that the reformulation of the problem has no influence on the convergence of the calculations.

If properly interpreted, the topology optimization results may be an inspiration for architects and engineers at an early design stage. One can see that for generated topologies, small holes or slender rods can occur. Those elements can be difficult to manufacture, and in the case of slender rods, they can be exposed to the loss of stability. To prevent the problems mentioned, either manufacturing constraints can be applied or those elements need to be modified or removed in the final design steps.

4.1.3. Topology Optimization of Multi-Material Structure under External Load and the Self-Weight

Next, the optimization process was performed for the same structure, but self-weight loading was additionally included, with the gravity acceleration equal to 9.81 m/s². The implementation of self-weight loading has a significant influence on the final results of the topology optimization process, which is very important, especially when dealing with the optimization of massive engineering structures, for example, bridges or carrying systems of tall buildings. Topology optimization of structures under design-dependent loads cannot be treated as an extension of standard formulation of topology optimization under fixed loads and therefore requires modifications of known and well-recognized approaches and algorithms. The presented study shows that the proposed local rule can be effective also for the more complicated problems like design-dependent loads. The final topology is shown in Figure 13. The value of compliance for the final topology equals 2.531 × 10⁻² Nm, and the final displacement equals 0.336 × 10⁻⁵ m, while for the initial structure they are 4.333 × 10⁻² Nm and 0.454 × 10⁻⁵ m, respectively.

For the above example (8005 elements), 30 iterations were needed to reach a convergence, which took about 3 min (one half of the structure was investigated due to the symmetry). The iteration history showing the convergence to the final solution is provided in Figure 14.

The calculations were repeated for the corresponding three-dimensional structure (thickness of 0.25 m), as was done in the previous section, and the final results are presented in Figure 15. The value of final compliance for the optimized structure equals 0.204 Nm (0.399 Nm for the initial structure), and the final maximal displacement equals 0.431 × 10⁻⁴ m (0.833 × 10⁻⁴ m for the initial structure).

For this example (20,808 elements), 35 iterations were needed to reach a convergence, which took less than 12 min (onequarter of the structure was investigated due to the symmetry). The iteration history showing the convergence to the final solution is provided in Figure 16.

As for the previous example, the considerations regarding the multi-material 2D structure were extended by replacing the volume constraint with the mass constraint. In what follows, the mass constraint µ = 0.5 was implemented for the same geometry, material data and boundary conditions. This time, since the multi-material structure was investigated, the constraint replacement influences the final results. The compliance of the optimized structure equals 0.276 × 10⁻² Nm, while the maximal displacement 0.138 × 10⁻⁵ m. The final topology is presented in Figure 17. The iteration history showing the convergence to the final solution is provided in Figure 18 (30 iterations were needed to get the solution). It is worth underlining that the reformulation of the problem has no influence on the convergence of the calculations.

4.2. Example 2—Multi-Layer Structure

A multi-layer structure presented in Figure 19 was chosen as the next example. This structure can model a massive engineering structure for which the consideration of self-weight loading is required. In addition, the application of multi-material solutions for this type of structure can be found in architectural and structural engineering practice. In this case, the threshold values c₁* = c₂* are selected and calculated based on the average compliance of the structure.

The multi-layer structure is built of two materials. Both materials are linear, elastic and isotropic. Material 1 (M₁ drawn in Figure 19 in dark grey color) is characterized by the Young modulus equal to E₁ = 2.1 × 10⁵ MPa, the Poisson ratio ν₁ = 0.3 and density ρ₁ = 7860 kg/m³; whereas material 2 (M₂ drawn in Figure 19 in light grey color) is characterized by the Young modulus equal to E₂ = 2.7 × 10⁴ MPa, Poisson ratio ν₂ = 0.2 and density ρ₂ = 2500 kg/m³. The black regions simulating girders are excluded from the design process. The structure is loaded by the external pressure equal to 2 10⁴ N/m acting on each girder as is shown in Figure 19.

4.2.1. Topology Optimization of Multi-Material Multi-Layer Structure under External Load and Self-Weight

The topology optimization of the structure under external was performed, and the results are presented in Figure 20a. For the final topology, the values of compliance and maximal displacement are 11.164 Nm and 0.242 × 10⁻⁴ m, while for the initial structure, they are equal to 17.791 Nm and 0.326 × 10⁻⁴ m, respectively.

Next, the generation of optimal topology for the multi-layer structure was performed, but this time with the self-weight loading included. The final topology is shown in Figure 20b. The value of compliance for the final topology equals 544.38 Nm, whereas the maximal displacement is equal to 0.191 × 10⁻³ m. For the initial structure, the respective values are 1379.33 Nm and 0.268 × 10⁻³ m.

For this example (15,300 elements), 30 iterations were needed to reach a convergence, which took about 4 min. Figure 21 shows the iteration history for the convergence to the final solution without self-weight loading (topology presented in Figure 20a). Figure 22 shows the solution with self-weight loading (topology presented in Figure 20b).

4.2.2. Topology Optimization of Multi-Material Multi-Layer Structure under External Load and Self-Weight—Design of the Reinforcement

The present consideration regards the implementation of the sub-domain-oriented topology optimization to the special case of the design domain selection. In what follows, only one part of the multi-layer structure, selected in three different ways, is considered to be the design domain. As previously, the external pressure and self-weight loading were applied. The final topologies obtained for the three cases under consideration are presented in Figure 23. One can observe that the location of the selected design domain significantly influences not only the generated topology but also the final value of compliance. The results of calculations are 2025.02 Nm for the structure on the left, 1756.25 Nm for the middle one and 1501.25 Nm for the one on the right. The above-mentioned approach allows for a sample modification of supports of existing structures with a view to obtain the minimal compliance of the reinforcement.

5. Topology Optimization Techniques Adapted to Strengthening and Retrofitting of Civil Structures

In the case of cultural and architectural heritage, strengthening weakened material and retrofitting structural damage can be a concern for historians and architects, but also for engineers. The implementation of innovations in engineering and science seems to be an attractive alternative to conventional, well-known ideas. Following this line of reasoning, this paper proposes a concept for optimal designing of a restoration process, which allows reducing the mass of strengthening elements, including one or more applied materials.

5.1. Strengthening of Structures Suffering from the Effects of Material Degradation

An attempt to build a simple but effective technique for design strengthening for cultural heritage monuments and historical buildings suffering from the effects of material degradation is made in this section. It is proposed to treat the restoration process as a topology-optimization task, in which the strengthening with a low volume of material and maximal stiffness are sought for.

As an example of strengthening civil structures suffering from the effects of material degradation using multi-domain and multi-material oriented topology optimization, the bridge structure presented in Figure 24 was selected. The bridge is loaded by a uniform pressure equal to 1000 N/m.

The structure is divided into two domains. The first one, a design-passive region, is presented in Figure 25, where grey color marks the weakened material. The design-active domain is marked red—in this domain, the strengthening is designed by topology optimization. For the first case, uniform material was selected for the whole structure (Young modulus E = 27 × 10⁹ Pa, Poisson ratio ν = 0.3 and density ρ = 2500 kg/m³), and then the randomly weakened material was selected for the design-passive region to model a material degradation in the damaged construction. The material’s random properties were calculated from the interval: 0.8 × 27 × 10⁹ Pa–27 × 10⁹ Pa and 0.8 × 2500 kg/m³–2500 kg/m³. These material uncertainties may serve as a simplified model of real-life material of old structures, for which it is impossible to describe material characteristics at every point without detailed examination.

The calculations for uniformly distributed material were performed first, which allowed for reducing the compliance value from 3.006 × 10⁻² Nm for the initial structure (without strengthening) to 1.871 × 10⁻² Nm for a structure with strengthening (see Figure 26a). At the same time, the maximum displacement was reduced from 0.315 × 10⁻⁵ m to 0.150 × 10⁻⁵ m. The calculations for random weakness of the material were performed next, and as the result, the initial value of compliance of the structure, 3.334 × 10⁻² Nm, was reduced to 2.061 × 10⁻² Nm (see Figure 26b). At the same time, the maximum displacement was reduced from 0.348 × 10⁻⁵ m to 0.164 × 10⁻⁵ m. What is worth underlining is that the volumes of the strengthening for both cases are the same (volume fraction κ for both cases equals 0.25).

As it was shown preliminarily in [49], the proposed approach can be an effective tool for designing the strengthening of damaged structures in the case of weakened material. In the present study, the self-weight loading has been taken into account, and the illustration of the influence of the self-weight on the final topologies is presented in Figure 27.

First, the calculations including self-weight for uniformly distributed material were performed, which allowed the compliance value to be reduced from 162 Nm for the initial structure to 109 Nm for a structure with strengthening (see Figure 27a). At the same time, the maximum displacement was reduced from 0.211 × 10⁻³ m to 0.105 × 10⁻³ m. Next, the calculations for a random weakness of the material including self-weight were performed, and the initial value of compliance of the structure 144 Nm was reduced to 100 Nm (see Figure 27b). At the same time, maximum displacement was reduced from 0.221 × 10⁻³ m to 0.104 × 10⁻³ m.

For example, as presented in Figure 26a, 60 iterations (about 5 min.) were needed to reach a convergence, while for the example presented in Figure 26b, 80 (about 7 min.). The iteration history showing the convergence to the final solution is provided in Figure 28 and Figure 29.

For examples where the self-weight loading is applied, 80 (i.e., about 7 min) and 50 (i.e., about 4 min) iterations were needed to reach the convergence. The iteration history showing the convergence to the final solution is provided in Figure 30 and Figure 31 (for topologies presented in Figure 27a,b, respectively).

An interesting observation is that the final topology for examples for which the self-weight was included is almost the same. The explanation can be found in the ratio of values of external loadings and self-weight: when the self-weight is applied and the external forces are low compared to it, the distribution of the material is governed by its own weight more than by external loads transmitted through the weakened material to the design area.

5.2. Strengthening of Structures Suffering from the Structural Damage—The Multiple-Load Case

As an example of strengthening civil structures suffering from the effects of structural damage, a model of a building facade is proposed (see Figure 32). The idea is to use the given material amount to build the strengthenings, which can minimize the compliance of the stucture. To illustrate the versatility of the method, a multiple-load case is investigated. The loading cases are presented in Figure 33. The loading models wind loads and a snow load. The horizontal pressure equals 600 N/m, and the vertical load is distributed, as is shown in Figure 33c, equal to 7000 N.

The applied material data are the Young modulus E = 2 × 10¹⁰ Pa, the Poisson ratio ν = 0.3. In these cases, the material data are the same for the structure and the strengthening. The volume fraction κ for all cases for design-active regions equals 0.2. Figure 34 presents the final topologies for all schemes of loading treated separately and, at the end, the final topology for a multiple-load case. The magnification of the design-active region is applied.

The final values of compliances are, for the final topology for load case I, 6.133 × 10⁻⁴ Nm; for the final topology for load case II, 6.262 × 10⁻⁴ Nm; for the final topology for load case III: 3.978 × 10⁻³ Nm. Considering the multiple-load case, the compliance should be checked for final topology, but for every loading separately, and what follows, the compliance values for the considered load cases are 0.692 × 10⁻³ Nm, 0.676 × 10⁻³ Nm and 4.032 × 10⁻² Nm, respectively.

For this example, the convergence of the multiple-load case presented in Figure 34d is discussed. At every iteration, three analyses are performed, so the computational time is larger than for the single-load case. For the presented example, 30 iterations were needed for the solution, which took about 9 min for 14,087 elements used for discretization. The iteration history showing the convergence to the final solution is provided in Figure 35.

5.3. Strengthening of the Multi-Material 3D Structures

The strengthening of a multi-material 3D structure was selected to present the capabilities of the method in application to topology optimization for real three-dimensional engineering structures. The model of breach as a model of structural damage is presented in Figure 36, with the final strengthening (drawn in Figure 36 in red color). The material data for the strengthening was defined as E₁ = 2 × 10¹¹ Pa and ν₁ = 0.3, while the material data for the damaged structure (drawn in Figure 36 in grey color) was defined as E₂ = 2 × 10¹⁰ Pa and ν₂ = 0.3. The volume fraction was defined as 0.4. The considered structure was divided into bottom (1 m × 0.5 m × 0.1 m) and top (1 m × 0.5 m × 0.1 m) blocks, determined as the non-optimized subdomains, and a middle block (1 m × 0.5 m × 0.8 m), determined as an optimized sub-domain. The bottom plane of the structure was fixed, while the top-right corner was loaded by three concentrated forces equal to 1000 N each, as presented in Figure 36. The initial compliance of the structure 4.990 × 10⁻² Nm was reduced by the optimization process to the value of 2.511 × 10⁻² Nm.

For the 3D example, a quarter of a structure was considered with 20,230 elements. A total of 30 iterations were needed to get the solution, which took about 8 min. The iteration history showing the convergence to the final solution is provided in Figure 37.

6. Discussion and Conclusions

This paper presents the idea of sub-domain-based topology optimization of multi-material structures under both external and self-weight loading. The design domain is divided into sub-regions, for which materials with different properties are applied. This approach to topology optimization allows one to obtain minimal-compliance structures with less material compared to the initial design. As has been shown, introducing a bi-material structure leads to different results when compared to the problem formulated for uni-material structure. On the other hand, the distribution of material within a sub-domain allows for this concept to be applied to architectural design, while multi-material solutions are desirable in innovative sustainable buildings to make constructions more environmentally friendly and modern. Furthermore, designers can be more flexible during the design process, with the possibility of using standard building material and new alternative solutions like renewable materials. The main issue, which can be a challenge in designing massive structures using topology optimization, is the implementation of the design-dependent loading, namely self-weight loading. In many cases, when applying optimization algorithms, numerical instabilities occur. The chosen method allows manufacture-friendly designs to be obtained easily and in a short calculation time, even when dealing with materials of different physical properties. Additionally, the proposed method is utilized for the sub-domain optimal design of strengthening of structures suffering from structural damage and/or material degradation. Material uncertainties, design-dependent loading (self-weight) and a multiple-load case are considered to illustrate the versatility of the method and the possibility of defining more complicated problems. What is very important is that the concept topology optimization performed by Cellular Automata is a fast-converged technique, not suffering from the checkerboard and mesh dependency effect. This method can be successfully used in practical civil engineering problems.