A hybrid topology optimization method applied to reinforced concrete structures using polygonal finite elements

Um método híbrido de otimização topológica aplicado em estruturas de concreto armado utilizando elementos finitos poligonais

Abstract

This work introduces a new alternative to obtain strut-and-tie models using the hybrid topology optimization method, which is already proposed in the technical literature and is refined here to use polygonal finite elements and accelerate the solution of the material nonlinearity problem. In this method, concrete is approached as a continuum, using polygonal two-dimensional finite elements, and steel bars as truss elements, using one-dimensional finite elements with two nodes. For a closer representation of reality, an orthotropic constitutive model for concrete was implemented considering different compression and tensile stiffness values, which is one of the advantages of the model. Further, the hybrid method limits the final layout of steel bars, thereby generating better structures from a constructive point of view, while allowing greater freedom for the shape and concrete strut slope. However, this method is more complex, and it increases the computational cost, which was substantially minimized through the implementation of an algorithm. Results obtained for some domains were very close to the results of other methodologies; however, small differences were noted that may be relevant to the final result. Other domains showed results with greater differences, thereby significantly changing the final strut-and-tie model and presenting a new structural design alternative.

Keywords:
reinforced concrete; topology optimization; hybrid; orthotropic; strut-and-tie

Resumo

Esse trabalho apresenta uma nova alternativa para a obtenção de Modelos de Bielas e Tirantes para estruturas de Concreto Armado através do Método de Otimização Topológica Híbrida, já proposto na literatura técnica, sendo aqui refinado para utilizar elementos finitos poligonais e acelerar a solução do problema da não linearidade material. Nesse método, o concreto é aproximado como meio contínuo, utilizando elementos finitos bidimensionais poligonais, e as barras de aço como elementos de treliça, utilizando elementos finitos unidimensionais com dois nós. Para uma representação mais próxima da realidade, implementou-se um modelo constitutivo ortotrópico para o concreto, considerando rigidezes distintas para a compressão e tração, sendo uma das vantagens do modelo. O método híbrido também limita a disposição final das barras de aço, gerando estruturas melhores do ponto de vista construtivo, enquanto permite uma maior liberdade para a forma e a inclinação das bielas de concreto. Como desvantagem, o método é mais complexo, aumentando o custo computacional, que porém foi substancialmente reduzido através da implementação de um procedimento. Os resultados obtidos para alguns domínios se aproximaram bastante dos resultados de outras metodologias, mas pequenas diferenças são notadas que podem ser relevantes para o resultado final. Outros domínios apresentaram resultados com maiores diferenças mudando de forma significativa o modelo de bielas e tirantes final, apresentando, assim, uma nova alternativa de dimensionamento.

Palavras-chave:
concreto armado; otimização topológica; híbrida; ortotrópico; bielas e tirantes

1 INTRODUCTION

Reinforced concrete is widely used in civil engineering because of its characteristics such as good resistance and versatility in obtaining structures with different shapes and because the construction techniques are well known. Concrete and steel work together, i.e., the adhesion between these two materials enables the transfer of forces between them. Despite being a widely used material, the behavior of reinforced concrete is complex.

At the end of the 19th century, Wilhem Ritter and Emil Mörsch developed approaches to describe the behavior of reinforced concrete considering that a cracked beam of this material behaved similar to a truss of concrete struts and steel ties, thereby idealizing the path of forces in the structure. This model is known as the model of Mörsch's truss analogy or the strut-and-tie model (STM). Since then, several studies have been conducted to perfect the theory about this model. However, traditional STMs could easily describe only type B regions—regions in which Bernoulli's hypotheses apply—but they were unable to describe type D regions—regions located in discontinuity zones where the theory of Bernoulli does not provide good approximations with adequate accuracy [¹], [²]. Figure 1 shows these types of regions.

Figure 1
Type B and D regions according to Schlaichet al. [²].

Thus, in 1987, Schlaichet al. [²] developed a methodology for determining STMs with sufficient efficiency for type D regions that guarantee with reasonable security —through the lower plasticity theorem—the adequate performance of structures for the ultimate limit state of rupture. Several valid models with different reinforcement positions can be obtained with this methodology; therefore, the importance of a geometric steel-concrete configuration that provides the greatest possible strength with the least use of material is remarkable. Therefore, optimization of the arrangement of steel reinforcement in concrete is desired.

Optimization implies obtaining the most effective solution by providing appropriate criteria and restrictions. In structural engineering, structural optimization is divided into three types: dimensional, shape, and topology. Topology optimization, which is the focus of this work, modify how elements of a given structure are connected, thereby allowing an adoption of null sections for the elements, i.e., these elements are removed from the final structure [³]. The topological optimization method is used to identify strut-and-tie arrangements that will provide greater stiffness to a reinforced concrete structure. Several studies have been conducted to obtain several STMs. Kumar [⁴] applied optimization concepts existing at the time to obtain these models; subsequently, several other studies were conducted in this field of research [¹], [⁵]–[⁹], which demonstrates their importance in structural engineering with a significant highlight of topology optimization. In principle, it consists of two methods commonly used for the determination of these models: topological optimization method of discrete medium (TOM-D) and the topology optimization method of continuous medium (TOM-C).

In the TOM-D model—known as the ground structure optimization method—the domain to be optimized is discretized as a dense mesh of bars interconnected by nodes as a truss structure, and this truss structure is optimized. The advantage of this method is the possibility of obtaining optimal structures with better conditions of applicability in practice. However, the disadvantage is that the loading will follow straight paths according to a pre-defined truss mesh, which will also cause a dependence on the adopted mesh as to the result.

The TOM-C model optimizes the continuous environment of the domain, and it uses two-dimensional or three-dimensional finite elements for its discretization and admits only the binary character with or without the presence of the material. The advantage of this method is the freer nature of the continuous structure that allows the path of forces to assume any shape and allows the members to connect at any angle. The disadvantage of this methodology is that this greater freedom can generate optimal structures that are impractical to build.

Although these methods have a good efficiency, new topology optimization methods with a specific focus on reinforced concrete have been developed; one such approach is proposed by Amir [⁸]. In their study, concrete is represented as an improved gradient model of continuous damage with stress relaxation, and the reinforcement is modeled by elastic bars inserted into the concrete domain using a ground structure. Therefore, steel bars and concrete are modeled differently as they behave differently in reinforced concrete. Following this, Gaynoret al. [¹] modeled concrete and reinforcement considering this distinct behavior. In their study, the continuous medium of the concrete is approximated by two-dimensional orthotropic finite elements using a four-node quadrilateral element; the steel is approximated by one-dimensional elements. The method adopted by Gaynoret al. [¹] was used in the work developed here, and it was modified to use polygonal isoparametric elements, which is referred to as the hybrid topology optimization method (TOM-H).

In this method, while the steel bars are approached by the truss bars, which are responsible for only tensile strength, the concrete is approximated by two-dimensional continuous finite elements, resisting compression stresses. Constitutive relationships used for concrete are obtained from the orthotropic model presented by Darwin and Pecknold [¹⁰]. This method has the following advantages: obtaining a more realistic topology for reinforced concrete structures than isolated methodologies (TOM-C and TOM-D) as the method considers the interaction of the two existing materials in reinforced concrete (i.e., concrete and steel) and it adopts an orthotropic constitutive model for concrete, which better represents its real behavior; providing greater freedom and better visualization of the flow of the compressive stress fields that represent the compressed concrete struts; and presenting a better practicality regarding the disposition of the steel. As a disadvantage, there is an increase in the computational cost in relation to the continuous and discrete isolated formulations; however, considering the current computational processors does not become a limiting factor to their use. Even so, in this work, a calculation procedure was elaborated that reduced this computational time considerably.

As a starting point, we used three algorithms written in MATLAB already existing in the technical literature by implementing the suggested methodology. These are:

• PolyMesher [¹¹]: Algorithm for creating meshes of convex polygonal two-dimensional elements using centroidal Voronoi diagram concepts.

• PolyTop [¹²]: Algorithm that calculates the optimal topological structure based on a given domain using two-dimensional polygonal elements, considering a fixed volume restriction and a formulation chosen by the user to penalize the density of the elements. This program was used as a basis for the development of this work.

• GRAND [¹³]: Algorithm that generates a truss mesh (ground structure) from a given domain; the volume of the truss mesh is optimized by the plastic formulation.

The algorithms were implemented in Python using the libraries Numpy [¹⁴], SciPy [¹⁵] Numba [¹⁶] and Matplotlib [¹⁷]. These algorithms inherited the name of the MATLAB version by adding the suffix py at the end to differentiate them. In addition, a new algorithm called HybridToppy [¹⁸] was developed based on the aforementioned algorithms and the hybrid method [¹], which solves optimization problems of two-dimensional reinforced concrete structures using convex polygonal finite elements and the hybrid methodology method.

2 PROPOSED FORMULATION

Optimization implies obtaining the best possible solution for a given problem while respecting the constraints imposed; mathematically, it is simply finding the minimum or maximum of a given function submitted to a series of constraints. In this study, we aim to find an optimal topology that provides the greatest possible stiffness of the structure that can be related to the total strain energy (compliance) defined by Equation 1.

$C\left(\mathit{\rho }\right)={\mathit{F}}^T\mathit{u}\left(\mathit{\rho }\right)$ (1)

Where $\rho$ denotes design variables: the area of the truss element and the density of the plate element for a unitary thickness; $\mathit{F}$ denotes the vector of the external forces applied to the structure; and $\mathit{u}\left(\mathit{\rho }\right)$ represents the vector with the nodal displacements of the structure.

The less the total strain energy, the greater is the stiffness of the structure; therefore, the problem is to minimize Equation 1 [³].

Displacements of the structure $u$ need to be determined to perform the optimization, and therefore, it is necessary to establish how the stress distribution in the structure is modeled. The finite element method (FEM) is generally used because it provides satisfactory results and has adequate compatibility with the existing optimization methods [¹], [⁵]–[⁹].

2.1 FEM formulation for steel reinforcement

The reinforcement in reinforced concrete was modeled as a one-dimensional axial element used in truss problems. Disregarding the load distributed in the reinforcement, it is possible to obtain the following formulation for the finite element corresponding to a bar with length $L$ in the local coordinates through FEM, expressed by Equation 2.

${\mathit{K}}_{\mathit{e}}^{\text{l}}{\mathit{u}}_{\mathit{e}}^{\text{l}}={\mathit{F}}_{\mathit{e}}^{\text{l}}$ (2)

Where (Equation 3)

${\mathit{K}}_{\mathit{e}}^{\text{l}}=\frac{ES}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$ (3)

is the local stiffness matrix of the element; ${\mathit{u}}_{\mathit{e}}^{\mathit{l}}={\left[\begin{array}{cc}{u}_{e1}^l& u_{e2}^l\end{array}\right]}^T$ denotes the local nodal displacement vector of the element; ${\mathit{F}}_{\mathit{e}}^{\mathit{l}}={\left[\begin{array}{cc}{F}_{e1}^l& F_{e2}^l\end{array}\right]}^t$ denotes the vector of the local nodal forces of the element; E denotes the elasticity modulus of the element; S denotes the cross-sectional area of the element; and L denotes the length of the element

Using a rotation matrix, it is possible to transform this system into the global coordinates of the structure. Then, a global stiffness matrix of the structure is created that encompasses all equations and degrees of freedom for each element. This matrix is produced by expanding the various local matrices, and by assigning zeros in degrees of freedom that are not part of the element. Thus, the equilibrium system is given by Equation 4.

$\mathit{K}\mathit{u}=\mathit{F}$ (4)

where $\mathit{K}$ = the global stiffness matrix of the structure, $\mathit{u}$ = the nodal displacement vector of the structure, and $\mathit{F}$ = the nodal force vector of the structure.

2.2 FEM Formulation for concrete

The concrete present in the structure is approximated as a domain discretized into two-dimensional finite elements. The problem was assumed to be flat, wherein the load belongs to the plane of the structure and there is no rigidity in the orthogonal direction. Thus, there are only two degrees of freedom per node, horizontal displacement, and vertical displacement. The formulation for the two-dimensional finite elements follows the same logic adopted for the one-dimensional elements.

The material constitutive law for this type of problem is given by Equation 5.

$\mathit{\sigma }=\mathit{D}\mathit{\epsilon }$ (5)

where

$\mathit{\sigma }={\left[\begin{array}{ccc}{\sigma }_x& {\sigma }_y& {\tau }_{xy}\end{array}\right]}^T$ and $\mathit{\epsilon }={\left[\begin{array}{ccc}{\epsilon }_x& {\epsilon }_y& {\gamma }_{xy}\end{array}\right]}^T={\left[\begin{array}{ccc}\frac{\partial u}{\partial x}& \frac{\partial v}{\partial y}& \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\end{array}\right]}^T$;

The constitutive matrix $\mathit{D}$ used in this work was obtained through Darwin and Pecknold’s [¹⁰] orthotropic model. In their work, the different behavior of concrete subjected to multiaxial stress is highlighted, and the proposed constitutive law defined in the main coordinates is given by Equation 6.

${\mathit{\sigma }}_{\mathit{p}}={\mathit{D}}_{\mathit{p}}\left({\mathit{\sigma }}_{\mathit{p}}\right){\mathit{\epsilon }}_{\mathit{p}}$ (6)

where ${\sigma }_p$ = the stress vector in the main coordinates; ${\mathit{D}}_{\mathit{p}}$ = the constitutive matrix in the main coordinates dependent on the main stresses; and ${\mathit{\epsilon }}_{\mathit{p}}$ = the strain in the main coordinates.

The constitutive matrix ${\mathit{D}}_{\mathit{p}}$ is defined by Equation 7.

${\mathit{D}}_{\mathit{p}}=\frac{1}{1-{\upsilon }_{eff}^2}\left[\begin{array}{ccc}{E}_1& {\upsilon }_{eff}{E}_{12}& 0\\ {\upsilon }_{eff}{E}_{12}& E_2& 0\\ 0& 0& \mathrm{0,25}\left(E_1+E_2-2{\upsilon }_{eff}{E}_{12}\right)\end{array}\right]$ (7)

Where the variables of matrix D_p are given by Equation 8.

$E_i=E_{ct}$, ${\upsilon }_i={\upsilon }_{ct}$ if ${\sigma }_{ci}>0$; $E_i=E_{cc}$, ${\upsilon }_i={\upsilon }_{cc}$ if ${\sigma }_{ci}\le 0$;

$E_{12}=\sqrt{E_1{E}_2}$; ${\upsilon }_{eff}=\sqrt{{\upsilon }_1{\upsilon }_2}$; ${\upsilon }_{ct}={\upsilon }_{cc}\frac{E_{ct}}{E_{cc}}$.(8)

Variables $E_{cc}$ and $E_{ct}$ are the elastic moduli of the concrete in compression and tensile, respectively. The variable ${\upsilon }_{cc}$ is the “Poisson's ratio of concrete in compression” adopted by Gaynoret al. [¹] as 0.2, and the same value is adopted in this study. The variable ${\upsilon }_{ct}$ is interpreted as the “Poisson's ratio of concrete in tensile.” The variable ${\upsilon }_{eff}$ is the effective Poisson's ratio that considers the stress acting in both directions.

Because the constitutive matrix is defined in the coordinates of the main stresses, it is necessary to calculate them and their respective directions. For this purpose, Equation 9 is used.

${\sigma }_{p\mathrm{1,2}}=\frac{{\sigma }_x+{\sigma }_y}{2}\pm \sqrt{{\left(\frac{{\sigma }_x-{\sigma }_y}{2}\right)}^2+{\tau }_{xy}^2}{\theta }_p=\frac{1}{2}{tan}^{-1}\left(\frac{2{\tau }_{xy}}{{\sigma }_x-{\sigma }_y}\right)$;(9)

Thus, it is possible to rotate Equation 6 from the main coordinates to the global coordinates, which is defined by Equation 5, using a rotation matrix in material coordinates. This rotation matrix can be deduced by rotating tensors and successive transformations to material coordinates. Therefore, the constitutive matrix for concrete is expressed by Equation 10.

$\mathit{D}={\mathit{Q}}^{\mathit{T}}{\mathit{D}}_{\mathit{p}}\mathit{Q}$ (10)

Where Q is given by Equation 11.

$\mathit{Q}=\left[\begin{array}{ccc}{cos}^2\left({\theta }_p\right)& {sen}^2\left({\theta }_p\right)& sen\left({\theta }_p\right)cos\left({\theta }_p\right)\\ {sen}^2\left({\theta }_p\right)& {cos}^2\left({\theta }_p\right)& -sen\left({\theta }_p\right)cos\left({\theta }_p\right)\\ -2sen\left({\theta }_p\right)cos\left({\theta }_p\right)& 2sen\left({\theta }_p\right)cos\left({\theta }_p\right)& {cos}^2\left({\theta }_p\right)-{sen}^2\left({\theta }_p\right)\end{array}\right]$ (11)

With the established constitutive relationship, the formulation of the FEM can be developed from the virtual work principle (VWP), which states that internal virtual work must be equal to external virtual work. When a body in balance is subjected to a field of virtual displacements that causes virtual deformations, and by using the appropriate matrix operations, the VWP can be described by Equation 12.

${\int }_{\Omega }\mathit{\delta }{\mathit{\epsilon }}^{\mathit{T}}\mathit{\sigma }d\Omega ={\int }_{\Omega }\mathit{\delta }{\mathit{u}}^{\mathit{T}}\mathit{f}d\Omega +{\int }_{S_{\sigma }}\mathit{\delta }{\mathit{u}}^{\mathit{T}}\mathit{q}dS$ (12)

where $\delta$ is the variational operator; $\mathit{u}={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$, where $u_1$ denotes the horizontal displacement and $u_2$ denotes the vertical displacement; $f$ refers to body forces; and $q$ refers to surface forces.

Based on this equation, it is possible to obtain a formulation that results in a relationship equal to that represented by Equation 2. Therefore, the stiffness matrix ${\mathit{K}}_{\mathit{e}}$ of the element and the nodal force vector ${\mathit{F}}_{\mathit{e}}$ of the element are given by Equation 13.

${\mathit{K}}_{\mathit{e}}=h_e{\int }_{{\Omega }_e}{\left(\mathit{B}\mathit{\varphi }\right)}^{\mathit{T}}\mathit{D}\left(\mathit{B}\mathit{\varphi }\right)d{\Omega }_e$ and ${\mathit{F}}_{\mathit{e}}=h_e{\int }_{{\Omega }_e}{\mathit{\varphi }}^{\mathit{T}}\mathit{f}d{\Omega }_e+{\int }_{S_{\sigma e}}{\mathit{\varphi }}^{\mathit{T}}\mathit{q}d{\Omega }_e$ (13)

Where $h_e$, ${\Omega }_e$, $\mathit{\varphi }$, and $\mathit{B}$ respectively denote the thickness of the element, domain represented by the area of the element, approximation function in the adopted finite element, and an operational matrix expressed by Equation 14.

$\mathit{B}={\left[\begin{array}{ccc}\frac{\partial }{\partial x}& 0& \frac{\partial }{\partial y}\\ 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]}^T$ (14)

The next step is to determine the type of two-dimensional finite element used for the approximation. In this work, the isoparametric polygonal finite element was used, which provides great flexibility for the generation of mesh, and it is better indicated for applications in the mechanics of solids that involve a significant change in topology in the domain material. Further, it presents better precision in numerical solutions and has less sensitivity to locking [¹⁹].

Talischiet al. [¹²] used the shape functions of Wachspress [²⁰], which is one of the options for shape functions of polygonal elements presented by Sukumar and Tabarraei [¹⁹]. These shape functions have desirable characteristics for these types of functions: being non-negative, interpolating, and forming a unit partition; having linear precision; and being precisely linear in the contour [¹⁹], [²¹]. To determine these shape functions and thus develop the finite element, it is interesting to use a regular polygon shape as a base. Interpolants can then be used both for approximating the displacement field through nodal values, and for transforming the coordinates of the initial element to the reference element of the integration. Figure 2 represents an example of an irregular hexagonal convex polygon and the related regular polygon showing points that facilitate the understanding of the process described below.

Figure 2
Polygonal element: (a) Irregular convex hexagonal polygon and displacement of its nodes. (b) Regular convex reference polygon.

The development of the Wachspress function is based on barycentric coordinates. Talischiet al. [¹²] used Sukumar and Tabarraei’s [¹⁹] definitions to establish the Equation 15 used in their work.

${\varphi }_i\left(\mathit{\xi }\right)=\frac{{\alpha }_i\left(\mathit{\xi }\right)}{{\sum }_{j=1}^n{\alpha }_j\left(\mathit{\xi }\right)}$ (15)

where ${\alpha }_i$ is an interpolating function expressed by Equation 16.

${\alpha }_i\left(\mathit{\xi }\right)=\frac{A\left(p_{i-1},p_i,p_{i+1}\right)}{A\left(p_{i-1},p_i,\mathit{\xi }\right)A\left(p_i,p_{i+1},\mathit{\xi }\right)}$ (16)

where $A$ represents the function that provides the area with the respective sign of the triangle formed by the points of the argument;$p$denotes the points referring to consecutive nodes of the polygon; and $\xi$ denotes the analyzed point and the approximation argument.

In a regular polygon, the parcel $A\left(p_{i-1},p_i,p_{i+1}\right)$ is the same for all $i,$ and it can be factored. Therefore, the notation of Equation 17 is adopted

$A_i\left(\mathit{\xi }\right):=A\left(p_{i-1},p_i,\mathit{\xi }\right)$ (17)

and Equation 18 is obtained.

${\alpha }_i\left(\mathit{\xi }\right)=\frac{1}{A_i\left(\mathit{\xi }\right)A_{i+1}\left(\mathit{\xi }\right)}$ (18)

The area can be given by the well-known formula of the matrix determinant for the area of a triangle, expressed by Equation 19.

$A_i\left(\mathit{\xi }\right)=\frac{1}{2}\left|\begin{array}{ccc}{\xi }_1& {\xi }_2& 1\\ p_{1,i-1}& p_{2,i-1}& 1\\ p_{1,i}& p_{2,i}& 1\end{array}\right|$ (19)

The partial derivatives of the area with respect to the reference coordinates are provided by Equations 20 and 21.

$\frac{\partial A_i}{\partial {\xi }_1}=\frac{1}{2}\left(p_{2,i-1}-p_{2,i}\right)$,(20)

$\frac{\partial A_i}{\partial {\xi }_2}=\frac{1}{2}\left(p_{1,i-1}-p_{1,i}\right)$,(21)

The derivatives of the interpolating functions are obtained by Equation 22.

$\frac{\partial {\alpha }_i}{\partial {\xi }_k}=-{\alpha }_i\left(\frac{1}{A_i}\frac{\partial A_i}{\partial {\xi }_k}+\frac{1}{A_{i+1}}\frac{\partial A_{i+1}}{\partial {\xi }_k}\right)$ (22)

Therefore, the derivatives of the shape functions with respect to the reference coordinates are given by Equation 23.

$\frac{\partial {\varphi }_i}{\partial {\xi }_k}=\frac{1}{{\sum }_{j=1}^n{\alpha }_j}\left(\frac{\partial {\alpha }_i}{\partial {\xi }_k}-{\varphi }_i{\sum }_{j=1}^n\frac{\partial {\alpha }_j}{\partial {\xi }_k}\right)$ (23)

By applying this function in Equation 13 and through numerical integration using a quadrature, the regular polygon is divided into triangles, each one solved individually and then summed to obtain the total integral. This process establishes a relationship equivalent to Equation 2. Adding all extended stiffness matrices, considering all degrees of freedom of the element, results in a relationship equivalent to that presented in Equation 4, which refers to the entire domain.

2.3 TOM-H Formulation

TOM-H has some peculiarities regarding the formulation of the FEM. Because reinforced concrete, concrete, and steel work together considering the same deformations and displacements to resist the imposed load, the stiffness of the structure is the result of the sum of the stiffnesses of the two materials, i.e., the stiffness matrix of the reinforced concrete is the sum of the stiffness matrix of the ground structure that represents steel, and the stiffness matrix of the continuous elements, that represents the concrete. Therefore, the meshes of the two types of elements must be compatible and have the same degrees of freedom. However, it is of interest that the reinforcing elements have less freedom of form than the concrete, thereby improving the constructive aspect; thus, the base mesh for the formation of the ground structure must be more sparse than the mesh of the two-dimensional continuous elements. Figure 3 represents the interaction of these two meshes, where it is observed that the base mesh used for the ground structure is much sparser than the mesh of continuous elements.

Figure 3
Mesh for solving the problem: (a) Base mesh for the generation Ground Structure; (b) Mesh of continuous elements; (c) Generated ground structure.

For the optimization in question, a bilinear constitutive model was considered for both materials, steel and concrete. For the steel material, only tensile stiffness was considered, because it is intended that the reinforcement acts only as an element of tensile strength. For the concrete material, in addition to compressive stiffness, a negligible tensile stiffness was considered to avoid singularities in the finite element; numerically, there cannot be degrees of freedom without the associated stiffness. The stress–strain relationship of the two materials is shown in Figure 4.

Figure 4
Graph of stress–strain relationship adopted in the hybrid optimization model. The stress measurement unit is irrelevant because what matters for the final STM is the ratio between elasticity modules.

The adoption of the bilinear constitutive model brings about a substantial improvement in representativeness, taking into account the distinct behavior of the two materials employed. Moreover, despite the nonlinearity of the final system generated by the finite element model, the model itself exhibits linear characteristics. This linear behavior not only simplifies the modeling process but also reduces the computational requirements, resulting in lower processing time and resource usage. As a result, there was no necessity to utilize a nonlinear finite element formulation. Furthermore, for design purposes, the bilinear constitutive model can be sufficient for the optimization process and the generation of a strut-and-tie model, eliminating the need for additional modeling complexities. Additionally, the methodology's outcome is more generic, preserving the principle of proportionality as long as the sign of the applied loads remains unaltered.

The optimization problem $\left(P\right)$ can be expressed by Equation 24.

$\left(P\right)=\left\{\begin{array}{c}{min}_{\rho }{\mathit{F}}^{T}u\left({\mathit{\rho }}_{\mathit{c}},{\mathit{\rho }}_{\mathit{t}},{\mathit{\sigma }}_{\mathit{c}},{\mathit{\sigma }}_{\mathit{t}}\right)\\ \text{s}\text{u}\text{b}\text{j}\text{e}\text{c}\text{t}\text{t}\text{o}{\mathit{\rho }}_{\mathit{t}}{\mathit{L}}_{\mathit{t}}+{\mathit{\rho }}_{\mathit{c}}{\mathit{A}}_{\mathit{c}}-V_{max}\le 0\\ {\mathit{\rho }}_{\mathit{t}}>\mathrm{0,0}\le {\mathit{\rho }}_{\mathit{c}}\le 1\end{array}\right.$ (24)

Where ${\mathit{\rho }}_{\mathit{c}},{\mathit{\rho }}_{\mathit{t}},{\mathit{A}}_{\mathit{c}},{\mathit{L}}_{\mathit{t}},V_{max}$ denote the density of the continuous elements, where zero indicates the absence of the material and one indicates a presence; the cross-sectional area of the truss bars, called the density of the bars; the area of the continuous elements; the length of the truss bars and the volume restriction adopted, respectively.

The problem is solved for a domain with unitary thickness. The modulus of elasticity of the steel and constitutive matrix depend on stress $\mathit{\sigma }$ in the structure; therefore, vector $u$ is implicitly provided by the nonlinear system $\mathit{K}\left(\mathit{\sigma }\right)\mathit{u}=\mathit{F}$ where the stiffness matrix $\mathit{K}\left(\mathit{\sigma }\right)$ is given by Equation 25.

$\mathit{K}\left(\mathit{\sigma }\right)={\mathit{K}}_{\mathit{t}}\left({\mathit{\sigma }}_{\mathit{t}}\right)+{\mathit{K}}_{\mathit{c}}\left({\mathit{\sigma }}_{\mathit{c}}\right)$ (25)

${\mathit{\sigma }}_{\mathit{t}},{\mathit{\sigma }}_{\mathit{c}}$ denote the stresses acting on the reinforcement and the concrete, respectively, and they are dependent on $\mathit{u}$; ${\mathit{K}}_{\mathit{t}}$ is obtained by the formulation of finite elements of truss; ${\mathit{K}}_{\mathit{c}}$ is obtained by the formulation of two-dimensional finite elements with a small change in the element stiffness matrix.

To adapt the formulation to topology optimization, which has a more binary view of the design variable, the presence or absence of material, ${\mathit{K}}_{\mathit{e}\mathit{c}}$(continuous element stiffness matrix) must be defined by Equation 26.

${\mathit{K}}_{\mathit{e}\mathit{c}}=\left(ϵ+\left(1-ϵ\right){\rho }_{ec}^q\right){\int }_{\Omega }{\left(\mathit{B}\mathit{\varphi }\right)}^T\mathit{D}\left(\mathit{B}\mathit{\varphi }\right)d\Omega$ (26)

where $ϵ\approx 0$ denotes a term to avoid singularity and $q$ denotes an exponent greater than or equal to one. Because ${\rho }_{ec}$ is only defined between zero and one, the intermediate values of the design variable will be penalized, thereby forcing the most binary view to be desired. However, this problem of pure topology optimization is ill-posed because if no restriction is established, it is always be possible to find an artificial topology even without physical applications such as a domain filled with bars as thin as we want, which will provide greater stiffness for the final structure [³]. The use of finite elements that imposes a restriction regarding the size of the mesh elements will present another ill-posed problem known as “Checkerboard” that consists of alternating values of the design variable without its convergence. This implies that a microscopic configuration on the mesh scale of the design variables provide similar responses on a macroscopic scale.

Thus, there is a need to apply a regularization filter, which avoids the sudden alternation of densities in nearby elements in the FEM approach of the two-dimensional elements, and therefore, it solves the illness. HybridToppy uses the same topological optimization formulation as Talischiet al. [¹²] presented below. The density of the continuous ${\rho }_c$ is now be provided by a convolution with a smooth filter function $\phi$, and can be expressed by Equation 27.

${\rho }_c\left(\mathit{p}\right)={\int }_{\Omega }\phi \left(\mathit{p},\stackrel{\text{\_}}{\mathit{p}}\right)x_c\left(\stackrel{\text{\_}}{\mathit{p}}\right)d\stackrel{\text{\_}}{p}$ (27)

where $\Omega$ denotes the domain of the structure, $\mathit{p}$ denotes the vector referring to the position of the point analyzed in the structure, and ${\mathit{x}}_{\mathit{c}}$ is be the problem variable for the concrete elements; however, it no longer has the physical meaning of density.

$\phi$ is given by Equation 28.

$\phi \left(\mathit{p},\stackrel{\text{\_}}{\mathit{p}}\right)={\left[{\int }_{B_r\left(\mathit{p}\right)\cap \omega }1-\frac{\left|\mathit{p}-w\right|}{r}dw\right]}^{-1}max\left(1-\frac{\left|\mathit{p}-\stackrel{\text{\_}}{\mathit{p}}\right|}{r},0\right)$ (28)

where $r$ denotes the filter radius used and $B_r\left(p\right)$denotes the circular domain of radius $r$ around $p$.

With the filter in place, there remains a need for one more procedure to solve the problem $P$ as the penalty has made the problem no longer convex, and it is now necessary to solve it through an iterative process with convex approximations. By defining $\mathit{x}$ as the problem variable that encompasses ${\mathit{\rho }}_{\mathit{t}}$ and ${\mathit{x}}_{\mathit{c}}$, the convex approximation of the problem (P) in interaction j in $x^j$ given by Equation 29 is obtained.

$\left(P^j\right)=\left\{\begin{array}{c}{min}_{x}C\left(x^j\right)+{\sum }_{i=1}^n{\left.\frac{\partial C}{\partial x_i}\right|}_{x_i^j}-\left(x_i^j\right)\left[{\left(\frac{x_i}{x_i^j}\right)}^{-1}-1\right]\\ \text{r}\text{e}\text{s}\text{t}\text{r}\text{i}\text{c}\text{t}\text{t}\text{o}g\left(x^j\right)+{\sum }_{i=1}^n{\left.\frac{\partial g}{\partial x_i}\right|}_{x_i^j}\left(x_i-x_i^j\right)\le 0\\ x_t>\mathrm{0,0}\le x_c\le 1\end{array}\right.$ (29)

where

$g\left(x\right)={\mathit{\rho }}_{\mathit{t}}{\mathit{L}}_{\mathit{t}}+{\mathit{\rho }}_{\mathit{c}}\left({\mathit{x}}_{\mathit{c}}\right){\mathit{A}}_{\mathit{c}}-V_{max}$ (30)

The convex approximation allows the use of the Lagrangian function and duality that, associated with the Karush–Kuhn–Tucker (KKT) conditions, allows solving the optimization problem $\left(P^j\right)$ through a numerical method of solving nonlinear equations. Thus, the iterative process provides the solution to the initial optimization problem (P). Talischiet al. [¹²] provided details of this procedure.

2.4 Algorithm used

Although HybridToppy exists as an independent algorithm, its use is correlated to several other algorithms. The first stage is performed in PolyMesherPy, in which two independent meshes are created, as shown in Figure 3, one for the continuous elements and one more for the formation of the ground structure. With the use of HybridMesherpy, these two meshes are made compatible, and therefore, the nodes of the base mesh of the ground structure coincide with the nodes of the mesh of the continuous elements. Subsequently, the ground structure is generated through the GenerateGS routine of GRANDpy.

Once the optimization process has been initiated, the structural problem is analyzed by the FEM, for which the equilibrium system expressed as a non-linear system needs to be solved iteratively. An alternative approach to solve this system using the fixed point logic was summarized by Gaynor et al. [¹]. However, it was noted that the stop criterion presented by them sometimes generated convergence problems using polygonal elements in HybridToppy because of minimal changes around the stress value of elements that could be considered zero. This caused a continuous change in the elasticity module of the materials, thereby taking too long to reach the stop criterion or never reaching it. These elements are negligible for the final result because all stresses are zero, from a numerical point of view. Thus, in this work, a stopping criterion for the non-linear analysis of the structural problem using the FEM, based on Equation 31, was used.

$\left(\frac{{\sum }_{j=1}^n\left|u_j^i-u_j^{i-1}\right|}{n}\right)÷\left(\frac{{\sum }_{j=1}^n\left|u_j^i\right|}{n}\right)<tol$ (31)

where $n$ denotes the number of degrees of freedom of the structure, $tol$ denotes a tolerance, and $i$ refers to the iteration analyzed.

The steps of the HybridToppy algorithm can be summarized as follows: resolve the nonlinear finite element system, update the variables using the optimality criteria, and check the convergence criteria of the optimization. The steps for solving the nonlinear finite element system include the calculation of local matrices of continuous elements, assembly of the global stiffness matrix of the continuous elements, assembly of the global stiffness matrix of the discrete elements, sum of the two matrices, system resolution, calculation of the stresses acting on all elements, saving variable values, and checking the finite element convergence criteria. If the convergence criteria are not met at any level, the algorithm returns to the first step of the corresponding level with the updated variables.

Finally, a graphical routine is available to display the optimization results in a visual format. This routine includes the option to apply a filtering mechanism, allowing the exclusion of bars with a relative area below a user-defined cutoff value, which in this study was set at 0.0001. Adjusting the cutoff values does not require re-running the entire optimization process and solely affects the graphical representation. However, caution must be used at this stage so that the use of the filter does not deviate from the optimal structural behavior proposed by the numerical result. For this, the user must pay attention to the number of bars that were excluded and the residue caused by this exclusion displayed by the routine.

The routine mentioned in the previous paragraph should be seen as the initial step in obtaining a strut-and-tie model. If the model generated, by analyzing the graphic result, is statically stable and using the inferior plasticity theorem, it is possible to guarantee the adequate performance of the structure.

2.5 Nonlinear Finite Element Method process accelerator

The basic procedure of the program was presented in a simple manner. It was observed that there was a considerable increase in the processing time of the program because of material nonlinearity, primarily due to the need to recalculate the stiffness matrices of the continuous elements. Further, for a tolerance factor of 0.5%, the number of material non-linearity iterations is reduced to 1 after a certain optimization iteration for each penalty value. This occurs because after a certain iteration, the optimization does not change the structure configuration significantly at each step.

$K^i$ is defined as the stiffness matrix of continuous elements in the optimization iteration i, and can be expressed by Equation 32.

${\mathit{K}}^{\mathit{i}}={\sum }_{e=1}^n{\mathit{K}}_{\mathit{e}}^{\text{*}\mathit{i}}\left({\mathit{\sigma }}_{\mathit{e}}^{\mathit{i}}\right){\rho }_e^i$ (32)

where ${\mathit{K}}_{\mathit{e}}^{\text{*}\mathit{i}}$ denotes the unit stiffness matrix considering the density with the unit value of an element $e$ in iteration i; ${\mathit{\rho }}_{\mathit{e}}^{\mathit{i}}$ denotes the density per area of the element $e$ in iteration i; and the $\sigma$ stresses in the element are dependent on the displacements $\mathit{u}$, which are dependent on the density distribution. However, after a certain number of optimization iterations, the structure configuration does not change considerably between one iteration and the other, the nature of the stresses—tensile or compression—remains almost the same for a small number of steps, and the angle $\theta$, which refers to the main directions of stresses, is slightly changed. Therefore, for a small number of $p$ steps, the variation of ${\mathit{K}}_{\mathit{e}}^{\text{*}}$ is small, which implied in Equation 33.

${\mathit{K}}_{\mathit{e}}^{\text{*}\mathit{i}}\approx {\mathit{K}}_{\mathit{e}}^{\text{*}\left(\mathit{i}+\mathit{p}\right)}$ (33)

Therefore, $K_e^{\text{*}}$ matrices can be considered as a constant by optimization iterations. The assembly of the unit stiffness matrix, which is a costly computational process, can be performed only for each determined number of iterations, and the ${\mathit{K}}_{\mathit{e}}$ matrices can be updated only by density, as calculated in the previous optimization iteration. Equation 34 expresses the stiffness matrix of continuous elements in the optimization iteration i+p.

${\mathit{K}}^{\left(\mathit{i}+\mathit{p}\right)}:={\sum }_{e=1}^n{\mathit{K}}_{\mathit{e}}^{\text{*}\mathit{i}}\left({\mathit{\sigma }}_{\mathit{e}}^{\mathit{i}}\right){\rho }_e^{\left(i+p\right)}$ (34)

This method was implemented so that the user can manually select how many iterations are necessary for the structure to stabilize, and only after applying the procedure described above, and the number of iterations p for which the local stiffness matrices of the continuous elements will not be calculated. For the adopted values of 10 and 5, a significant decrease in computational time was observed without prejudice to convergence or deviation of results.

3 RESULTS AND DISCUSSIONS

Four examples from the literature as shown in Figure 5 were evaluated to validate the method and compare the results obtained with those presented in the literature (Figure 6) and those obtained using other methodologies. Thus, the PolyToppy and GRANDpy algorithms were used to solve optimization problems with the TOM-C and TOM-D, respectively to compare the results obtained with HybridToppy.

Figure 5
Evaluated examples from literature: (a) Simple supported beam domain; (b) L-shaped wall beam domain; (c) Bridge beam domain; and (d) Shear wall with holes domain.

Figure 6
Results of other literature: black represent the concrete struts and red represent the steel ties; (a) Result of the domain simple supported beam by Gaynoret al. [¹]; (b) Result of the domain L-shaped wall beam by Gaynoret al. [¹]; (c) Result of the domain L-shaped wall beam by Zhonget al. [⁹]; (d) Result of the domain bridge beam by Gaynoret al. [¹]; (e) STM of the domain shear wall abstracted from the flow of elastic stresses by the FEM, according to FIB [²²].

The algorithm takes as input the number of elements in each domain, which determines the number of nodes and bars in the structure. Additionally, the constraint volume is considered, which, despite being a user choice, is a value that generates a result displaying a topology that can be associated with an STM. These main parameters for each analyzed domain are listed in Table 1.

Table 1
- Main parameters for each analyzed domain

Domain	Nodes	Elements on the Thin Mesh	Elements on the Coarse Mesh	Initial Bars	Volume Fraction (%)
Simple supported beam with poligonal base mesh	5987	3000	32	1822	50
Simple supported beam with regular base mesh	5987	3000	32	1822	50
L-shaped wall beam	3993	2000	42	1686	30
Bridge beam	3978	2000	30	2483	30
Sheared wall with holes	3992	2000	81	2095	30

3.1 Domain simple supported beam

This domain consists of a double-based beam that is 4 m long and 1 m high with a load concentrated in the middle, for which only half was discretized because of its symmetry, thereby saving computational time. The resulting topologies are shown in Figure 7.

Figure 7
Domain simple supported beam: (a) Result of GRANDpy; (b) Result of PolyToppy; (c) Result of HybridToppy with regular basis mesh; and (d) Result of HybridToppy with polygonal base mesh.

The GRANDpy and Polytoppy programs display the results in the same way of the original Matlab programs. The result displayed by GRANDpy (Figure 7a) shows bars topologies that optimize the problem without differentiating whether the stress is compressive or tensile. Meanwhile, Polytoppy provides a result (Figure 7b) that demonstrates the two-dimensional topology that optimizes the problem considering an isotropic material and given the volume constraint, where black symbolizes presence of material.

The result obtained by HybridToppy (Figure 7c, d) is represented by an optimized topology composed of steel bars, with colors ranging from blue to red, and concrete represented by the color black both with the stress–strain relationship presented by the Figure 4. The color bar on the right side of the figure refers to the cross-sectional area of the bars normalized, where the value 1 represents the section with the greatest area.

There was a similarity between the topologies of the three algorithms wherein the load is transmitted to the supports through elements dispersed in the form of a fan with transversal elements. The nature of the stresses in these elements is observed directly in the two topologies obtained with HybridToppy, the regular base mesh (Figure 7c), and the polygonal base mesh (Figure 7d).

For this example, the use of two types of base meshes for the generation of the polygonal and regular ground structure with square elements but with the same average element size were evaluated. However, there was no significant influence on the result obtained according to the type of the base mesh adopted.

Further, HybridTop provided a topology that resembles that presented by Gaynor et al. [¹] (Figure 6a).

3.2 Domain L-shaped wall beam

This domain was extracted from Kuchmaet al. [²³] with some adjustments made to the dimensions to facilitate the discretization of the ground structure with a regular mesh. Figure 8 presents the topologies obtained with the three algorithms.

Figure 8
Domain L-shaped wall beam: (a) Result by GRANDpy; (b) Result by PolyToppy; and (c) Result by HybridToppy.

The topologies of the three algorithms follow a similar trajectory with slight differences, wherein the load is dispersed in the region above the void to be transmitted to the two lower supports. The reinforcement zones are present above the orifice and in the lower region of the structure, and the compressed elements are abstracted from HybridToppy. The result obtained by HybridToppy is similar to that of Gaynor et al. [1] (Figure 6b).

Zhonget al. [⁹] conducted a more detailed analysis of this domain. In their work, the importance of a primary analysis was emphasized, wherein the STMs must correspond with the result of the linear stress analysis. Further, a subsequent analysis, with simulations, in which the correspondence of the reinforcement with the crack region in concrete was evaluated. They concluded that there is a region on the left edge of the domain where these analyses indicate the need for reinforcement that is not presented by the hybrid methodology. Therefore, Zhonget al. [⁹] elaborated their own STMs as presented in Figure 6c.

A justification for this possible difference was predicted by Bruggi [⁶]. He pointed precautions that must be employed to ensure the validity of the STMs obtained by the optimization methods. In his work, it was observed that the STMs resulting from his methodology needed additional reinforcement that would provide sufficient ductility to allow the use of the lower limit of plasticity theory. The optimization models adopted here do not present a restriction on deformation, thereby theoretically allowing the appearance of regions with great deformation, and consequently, excessive cracking, which is located where elements of the resulting topologies have zero densities (voids). Such regions, if any, must be reinforced with steel to ensure reasonable security, proper functioning of the structure respecting the ductility criteria, and requirements of the limit state of use.

3.3 Domain bridge beam

The bridge beam domain is a generic domain that resembles cross sections used in bridges. As in the simple supported beam, it is possible to take advantage of the symmetry of the structure, thereby saving computational time. The results obtained are shown in Figure 9.

Figure 9
Domain bridge beam: (a) Result of GRANDpy; (b) Result of PolyToppy; and (c) Result by HybridToppy.

There is a need for steel reinforcement in the upper part of the structure. In the lower zone, the structure is subjected to compression only, which is intuitively expected. There is a great similarity between the results obtained in all methodologies; however, it is noted that GRANDpy provides a result that requires a larger treatment. The results of HybridToppy and PolyToppy were quite similar, and based on the simplification adopted, the STMs abstracted from the two methodologies can be practically same. The result obtained from HybridToppy was also considerably similar to that presented by Gaynoret al. [¹].

3.4 Sheared wall with holes domain

The shear wall with holes domain was inspired by a domain exposed in the FIB report [²²]. The results obtained are shown in Figure 10.

Figure 10
Domain shear wall with orifice: (a) Result by GRANDpy; (c) Result by PolyToppy; and (d) Result by PolyTrussToppy with regular mesh.

The results of GRANDpy, PolyToppy, and HybridToppy show certain similarities; however, unlike others, HybridToppy presented a compression strut above the last hole. The constitutive model adopted by HybridToppy is closer to reality, as it considers the presence of two different materials (steel and concrete) and concrete as an orthotropic material, which is probably a more representative result.

For the complexity of the results, it is possible to note that the one provided by PolyToppy is the simplest; therefore, the treatment for the generation of an STM will be easier. The result of GRANDpy has an intermediate complexity and that of HybridToppy apparently has greater complexity. Thus, although the HybridToppy has a formulation closer to reality, it generates a more complex result that will require further treatment result to generate a strut and tie model.

The result provided by FIB [²²] (Figure 6e) differs from the solutions presented here because it was prepared without the optimization processes following Schlaich Schäfer and Jennewein's original methodology for STMs. As a disadvantage, the result does not necessarily present a model with optimum performance; however, the advantage is that it was established such that it greatly facilitates the execution as it is a simpler solution from a practical point of view.

4 FINAL REMARKS AND CONCLUSION

TOM-H presents itself as an alternative to traditional TOM-C and TOM-D methods, allowing greater freedom for concrete struts while restricting the steel reinforcement's ties, thereby improving the constructive aspect of steel reinforcement.

In this study, an algorithm based on the TOM-H and polygonal finite elements for the optimization of flat reinforced concrete structures was implemented. To represent the behavior of this composite material, the algorithm considered two types of materials: steel and concrete, with different tensile and compression stiffnesses. Further, concrete was approximated as an orthotropic material, thereby leading to a more complex formulation that was closer to reality. Thus, the insertion of material nonlinearity resulted in a higher computational cost; however, it was substantially reduced with the implementation of a procedure that reduced the disparity between the computational times of the methods.

Although the TOM-H constitutive model better describes the behavior of reinforced concrete, most of the domains analyzed by this methodology present results similar to the traditional ones, this shows that the traditional ones are adequate approximations for most cases. But the result from the Sheared Wall With Holes demonstrates that there will be domains in which a hybrid methodology will present different results from the traditional ones, which therefore must be analyzed in order to obtain more adequate strut and tie models.

The presented method is very promising for the study of reinforced concrete, and it allows efficient and direct determination of both the regions and optimal directions of reinforcement and the compression regions that most influence the stiffness of the structure. Therefore, it allows developing STMs that are likely to provide greater structural stiffness. For future research, it is suggested to refine this method by considering ductility criteria and controlling the cracking of the structure.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the Brazilian government agency CAPES, for the MSc scholarship of first author (Financing Code 001), and the Federal University of Sergipe of Brazil for the support in the development of this research. The conclusions and opinions expressed in this paper are the authors’ responsibility and do not reflect the points of view of the sponsors.

• Financial support: This work was carried out with the support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) Brazil – Financing Code 001.

• Data Availability:

  The data that support the findings of this study are openly available in HybridToppy, at https://github.com/pauloxvcr/hybridtop, or are available from the corresponding author, RMFC, upon reasonable request.
• How to cite: P. V. C. Rodrigues, and R. M. F. Canha, "A hybrid topology optimization method applied to reinforced concrete structures using polygonal finite elements," Rev. IBRACON Estrut. Mater., vol. 17, no. 3, e17306, 2024, https://doi.org/10.1590/S1983-41952024000300006

REFERENCES

• ¹
  A. T. Gaynor, J. K. Guest, and C. D. Moen, "Reinforced concrete force visualization and design using bilinear truss-continuum topology optimization," J. Struct. Eng., vol. 139, no. 4, pp. 607–618, Apr 2013, http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0000692
  » http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0000692
• ²
  J. Schlaich, K. Schäfer, and M. Jennewein, "Towards a consistent design of structural concrete," PCI J., vol. 32, no. 3, pp. 74–150, May 1987.
• ³
  P. W. Christensen and A. Klarbring, An Introduction to Structural Optimization. Dordrecht: Springer, 2009.
• ⁴
  P. Kumar, "Optimal force transmission in reinforced concrete deep beams," Comput. Struct., vol. 8, no. 2, pp. 223–229, 1978.
• ⁵
  Q. Liang, Y. Xie, and G. Steven, "Topology optimization of strut-and-tie models in reinforced concrete structures using an evolutionary procedure," ACI Struct. J., vol. 97, no. 2, pp. 322–330, 2000, http://dx.doi.org/10.14359/863
  » http://dx.doi.org/10.14359/863
• ⁶
  M. Bruggi, "Generating strut-and-tie patterns for reinforced concrete structures using topology optimization," Comput. Struc., vol. 87, no. 23–24, pp. 1483–1495, Dec 2009, http://dx.doi.org/10.1016/j.compstruc.2009.06.003
  » http://dx.doi.org/10.1016/j.compstruc.2009.06.003
• ⁷
  D.-K. Lee, C.-J. Yang, and U. Starossek, "Topology design of optimizing material arrangements of beam-to-column connection frames with maximal stiffness," Sci. Iran., vol. 19, no. 4, pp. 1025–1032, Aug 2012., http://dx.doi.org/10.1016/j.scient.2012.06.004
  » http://dx.doi.org/10.1016/j.scient.2012.06.004
• ⁸
  O. Amir, "A topology optimization procedure for reinforced concrete structures," Comput. Struc., vol. 114–115, pp. 46–58, Jan 2013, http://dx.doi.org/10.1016/j.compstruc.2012.10.011
  » http://dx.doi.org/10.1016/j.compstruc.2012.10.011
• ⁹
  J. T. Zhong, L. Wang, P. Deng, and M. Zhou, "A new evaluation procedure for the strut-and-tie models of the disturbed regions of reinforced concrete structures," Eng. Struct., vol. 148, pp. 660–672, Oct 2017, http://dx.doi.org/10.1016/j.engstruct.2017.07.012
  » http://dx.doi.org/10.1016/j.engstruct.2017.07.012
• ¹⁰
  D. Darwin and D. A. Pecknold, "Nonlinear biaxial stress-strain law for concrete," J. Struct. Eng., no. Apr, pp. 229–241, 1977.
• ¹¹
  C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes, "PolyMesher: a general-purpose mesh generator for polygonal elements written in Matlab," Struct. Multidiscipl. Optim., vol. 45, no. 3, pp. 309–328, Mar 2012, http://dx.doi.org/10.1007/s00158-011-0706-z
  » http://dx.doi.org/10.1007/s00158-011-0706-z
• ¹²
  C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes, "PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes," Struct. Multidiscipl. Optim., vol. 45, no. 3, pp. 329–357, Mar 2012, http://dx.doi.org/10.1007/s00158-011-0696-x
  » http://dx.doi.org/10.1007/s00158-011-0696-x
• ¹³
  T. Zegard and G. H. Paulino, "GRAND — Ground structure based topology optimization for arbitrary 2D domains using MATLAB," Struct. Multidiscipl. Optim., vol. 50, no. 5, pp. 861–882, Nov 2014, http://dx.doi.org/10.1007/s00158-014-1085-z
  » http://dx.doi.org/10.1007/s00158-014-1085-z
• ¹⁴
  T. E. Oliphant, A guide to NumPy, vol. 1. USA: Trelgol Publishing, 2006.
• ¹⁵
  P. Virtanenet al., "SciPy 1.0: fundamental algorithms for scientific computing in Python," Nat. Methods, vol. 17, p. 261-272, 2020, http://dx.doi.org/10.1038/s41592-019-0686-2
  » http://dx.doi.org/10.1038/s41592-019-0686-2
• ¹⁶
  S. K. Lam, A. Pitrou, and S. Seibert, "Numba: a LLVM-based Python JIT compiler," in Proc SC15: The Inter. Conf. High Perf. Comp., Network., Stor. Anal., New York, NY, USA, 2015, pp. 1-6. http://dx.doi.org/10.1145/2833157.2833162
  » http://dx.doi.org/10.1145/2833157.2833162
• ¹⁷
  J. D. Hunter, "Matplotlib: a 2D graphics environment," Comput. Sci. Eng., vol. 9, no. 3, pp. 90–95, 2007, http://dx.doi.org/10.1109/MCSE.2007.55
  » http://dx.doi.org/10.1109/MCSE.2007.55
• ¹⁸
  P. V. C. Rodrigues and R. M. F. Canha, ‘HybridTop.’ GitHub, Sep. 2020 [Online]. Available: https://github.com/pauloxvcr/hybridtop (accessed Mar. 1, 2023).
  » https://github.com/pauloxvcr/hybridtop
• ¹⁹
  N. Sukumar and A. Tabarraei, "Conforming polygonal finite elements," Int. J. Numer. Methods Eng., vol. 61, no. 12, pp. 2045–2066, Nov 2004, http://dx.doi.org/10.1002/nme.1141
  » http://dx.doi.org/10.1002/nme.1141
• ²⁰
  E. L. Wachspress A rational finite element basis. New York: Academic Press, 1975.
• ²¹
  G. Dasgupta, "‘Interpolants within Convex Polygons: wachspress’ Shape Functions’," J. Aerosp. Eng., vol. 16, no. 1, pp. 1–8, Jan 2003, http://dx.doi.org/10.1061/(ASCE)0893-1321(2003)16:1(1)
  » http://dx.doi.org/10.1061/(ASCE)0893-1321(2003)16:1(1)
• ²²
  Fédération Internationale du Béton. Task Group 4.4, Practitioners’ Guide to Finite Element Modelling of Reinforced Concrete Structures Lausanne, Switzerland: FIB - Féd. Int. du Béton, 2008.
• ²³
  D. Kuchma, S. Yindeesuk, T. Nagle, J. Hart, and H. Lee, "Experimental validation of strut-and-tie method for complex regions," ACI Struct. J., vol. 105, pp. 578–589, Sep 2008.

Edited by

Editors: Bernardo Horowitz, Guilherme Aris Parsekian.

Data availability

The data that support the findings of this study are openly available in HybridToppy, at https://github.com/pauloxvcr/hybridtop, or are available from the corresponding author, RMFC, upon reasonable request.

Publication Dates

• Publication in this collection
  04 Sept 2023
• Date of issue
  2024

History

• Received
  04 Mar 2023
• Accepted
  31 May 2023