Optimization‐based three‐dimensional strut‐and‐tie model generation for reinforced concrete

Strut‐and‐tie modelling (STM) is an effective and widely used method to design disturbed regions (D‐regions) of reinforced concrete structures. Among the various steps of STM, finding a suitable truss‐analogy model is the most challenging part. Even for experienced engineers it is difficult to find representative models for complex D‐regions, and this task is even harder for three‐dimensional (3D) D‐regions. To date, only a few 3D STM models have been proposed by researchers for several complex D‐regions, which leaves practitioners with little guidance. In order to solve this problem, a method is proposed to automatically generate 3D optimization‐based STM (3D‐OPT‐STM) models. The generation method comprises a compliance‐based topology optimization process that generates optimized structural forms by maximizing stiffness, a topology extraction method, and a shape optimization method. In this study, three 3D‐OPT‐STM models are generated for typical 3D D‐regions, and their performances are compared to manually created STM models. The generated 3D‐OPT‐STM models result in more economical designs. Moreover, geometrical and loading parameter studies demonstrate the applicability and robustness of the proposed method.


INTRODUCTION
It is a challenge for engineers to design disturbed regions (D-regions) of reinforced concrete (RC) structures due to their highly nonlinear strain distributions. The strut-andtie modeling (STM) has been established as an effective tool to design D-regions. STM was first generalized by Schlaich, Schafer, and Jennewein (1987) and Schlaich and Schafer (1991) as a versatile tool for designing RC structures. In the STM analysis, the flow of forces is idealized by abstracting concrete struts, reinforcement ties, and their connecting nodes. Using truss-analogy models that indicate the force and stress distribution, STM improves However, in current codes design guidances are provided for relatively simple D-regions only, such as deep beams. For relatively complex D-regions, current codes can only give limited directions. This evidently applies even more to disturbed three-dimensional (3D) domains than to disturbed 2D domains.
In the STM method, the first step is to simplify D-regions by discrete truss-analogy models to represent the force transfer mechanism. Various models can be used to design a single D-region and they can lead to significantly different performances (El-Metwally & Chen, 2017). It has been broadly accepted that finding a suitable truss-analogy model is an important part to guarantee that STM models will perform well (Bruggi, 2009;Kuchma, Yindeesuk, Nagle, Hart, & Lee, 2008;Liang, Xie, & Steven, 2000;Park, Yindeesuk, Tjhin, & Kuchma, 2010;Schlaich & Schafer, 1991). Although all STM models result in conservative designs, with an ill-chosen truss-analogy model the resulting design can be overly conservative, that is, uneconomical in terms of its steel usage.
Several approaches have been applied in the past decades to find suitable STM models. For example, load path methods were used to determine suitable trussanalogy models in Palmisano and Elia (2015) and Schlaich and Schafer (1991). Using stress field approaches, Schlaich and Schafer (1991) suggested conducting a linear finite element analysis and forming STM models based on the calculated stress field. Muttoni, Ruiz, and Niketic (2015) applied a similar concept. With increasing complexity of D-regions, the difficulties of generating suitable STM models also increase. In this regard, the use of topology optimization (TO) methods for STM model generation has become one of the most popular directions in this field. Researchers have proposed various TO methods for STM, for example, the evolutionary structural optimization (ESO) method (Almeida, Simonetti, & Neto, 2013;Kwak & Noh, 2006;Liang et al., 2000), the isoline method (Victoria, Querin, & Marti, 2011), the fullhomogeneous method (Herranz, Maria, Gutiérrez, & Riddell, 2012), and the SIMP (solid isotropic material with penalization) method (Bruggi, 2009(Bruggi, , 2010(Bruggi, , 2016Jewett & Carstensen, 2019). Optimization techniques have also been widely applied in the structural design beyond STM. Various applications have been reported in the field of the civil engineering. To name a few examples, Vantyghem, De Corte, Shakour, and Amir (2020) applied a TO method together with 3D printing to produce post-tensioned concrete structures. Stromberg, Beghini, Baker, and Paulino (2012) used continuum-beam integrated TO method to design laterally braced frames. Aldwaik and Adeli (2014) used a genetic algorithm to design high-rise buildings. Talischi, Paulino, Pereira, and Menezes (2010) applied the polygonal element in the TO process. Shape optimization techniques have been applied in designing truss structures (Adeli & Balasubramanyam, 1987Su, Wu, Ji, & Shen, 2018) and shell structures (Kociecki & Adeli, 2014Rombouts, Lombaert, De Laet, & Schevenels, 2019;Xia, Wu, & Hendriks, 2019).
Up to date, most investigations of finding STM models have focused on 2D D-regions. Compared to 2D STM models, there are much fewer studies on 3D STM models. One main reason is the increased difficulty of determining suitable STM models in 3D. However, as a result of the increased complexity of 3D D-regions, finding suitable 3D STM models is even more important. Cai (2002) proposed a 3D indeterminate STM model for a footing system. Leu, Huang, Chen, and Liao (2006) applied a refined ESO TO method to obtain optimized material layouts for creating 3D STM models. Similarly, based on the ESO TO methods, He and Liu (2010) generated STM models for anchorage diaphragms, Shobeiri and Ahmadi-Nedushan (2017) applied a bidirectional ESO method to obtain optimized material distributions for several 3D D-regions, and Zhong, Wang, Zhou, and Li (2017) used a microtruss-based ESO method for 3D STM models. Yun, Kim, and Ramirez (2018) generated 3D STM models using a grid-based approach. Dey and Karthik (2019) analyzed several pile-cap designs based on 3D compatible STM models. Several investigations (Meléndez, Sagaseta, Sosa, & Rubio, 2019;Souza, Kuchma, Park, & Bittencourt, 2009;Yun, Chae, & Ramirez, 2019) were conducted to investigate suitable 3D STM models for pile cap design problems.
Among different approaches, TO-based methods have the potential to provide a systematic procedure to generate suitable STM models. However, in current state, TO approaches only provide optimized material distributions as inspirations for STM models. Manual and subjective interpretations and adjustments are required to transform TO results (material distributions) to STM models (trussanalogy models). On the one hand, these subjective interventions hamper the systematic generation of STM models; on the other hand, the created STM models may not exhibit the required axial force equilibrium state. These problems have been demonstrated in our previous review (Xia, Langelaar, & Hendriks, 2020b). Therefore, in order to fill this gap and move towards developing a systematic STM method, an automatic generation method to generate suitable truss-analogy models for 3D STM is proposed in this paper. To the best of our knowledge, it represents the first method of its kind to address 3D STM model generation.
The method is developed based on our previously proposed optimization-based STM (OPT-STM) generation method (Xia, Langelaar, & Hendriks, 2020a), which was focused on 2D problems only. The proposed method in this paper includes three main steps: (a) the 3D compliance-minimization based TO, (b) the topology extraction that automatically generates truss-analogy models based on 3D TO results, and (c) the shape optimization that generates valid truss-analogy models for STM. In this paper the framework of the method is repeated briefly for convenience, while highlighting the modifications made for making the method applicable, effective, and robust for 3D D-regions. Key innovations this paper introduces in comparison to the 2D method to enable the treatment of 3D cases are: • various measures reducing the computation costs for the computationally demanding 3D TO process; • proposing a novel robust procedure in extracting 3D truss-like structures; • extending the STM shape optimization procedure for 3D models, including a continuous minimum length constraint.
Through the integrated optimization process, 3D-OPT-STM models can be generated in an automated manner, which simplifies the design process in engineering practice. Thus a systematic method for generating suitable 3D STM models for 3D D-regions is devised. Based on the proposed generation method, 3D-OPT-STM models for three typical D-regions (four-pile cap, corbel, and box girder) are generated and compared with the manually created STM models. The proposed generation method enables us to further investigate three factors in the STM method systematically: • A geometrical parameter study reveals the transition from one optimal topology to another as 3D aspect ratios are varied. • A loading parameters study shows the effect of adding out-of-plane load components on optimal transfer mechanisms. • Discretizing distributed loads is an essential step in the STM method. The effect of different load discretization schemes of surface loads in 3D is investigated.
Along with the variation studies, the applicability and robustness of the 3D-OPT-STM is demonstrated. It is emphasized that practical requirements, such as steel detailing and constructability aspects, are outside the scope of the current study.
The remainder of the paper is organized as follows. The 3D-OPT-STM generation method is introduced in Section 2. 3D-OPT-STM models for three typical D-regions are generated in Section 3. The performance of the generated STM models is evaluated by nonlinear finite element analysis (NLFEA) simulations. In Section 4, the three factors of the 3D STM are investigated by adopting the generation method. Finally, the conclusion of this paper is given in Section 5.

AUTOMATIC GENERATION METHOD FOR THREE-DIMENSIONAL OPTIMIZATION-BASED STRUT-AND-TIE MODELS
In this section, a generation method is proposed to automatically generate 3D-OPT-STM models, which forms a generalization of our previously proposed method limited to 2D cases (Xia et al., 2020a). The proposed method comprises three steps: (a) compliance minimization based 3D TO, (b) topology extraction, and (c) STM-based shape optimization. The optimized material layouts for various design problems can be obtained in the TO phase. Next, in the topology extraction phase, the obtained 3D TO results are transformed to truss-like structures that comprise a network of connected straight bars. Because the extracted truss-like structures are usually unstable structures, 3D beam elements are used to analyze their force distributions. In order to satisfy the requirement of axial force equilibrium for STM analysis, 3D shape optimization is conducted for the truss-like structures to generate adequate STM models in the last step. The flowchart of the proposed method is presented in Figure 1. The extensions and improvements of each step are presented in the following subsections.

2.1
Step 1: 3D TO optimization TO provides numerical material distributions for conceptual designs without specific shapes. Various methods and applications of TO methods have been reported after the initial work by Bendsøe and Kikuchi (1988). In this paper, we apply the classical and popular SIMP (solid isotropic material with penalization) method (Bendsøe & Sigmund, 2003) for the compliance minimization problem. The adopted SIMP method for 3D problems is a natural extension of our previous 2D work (Xia et al., 2020a), where the SIMP parameter settings were introduced and their influences on the generated STM models have been discussed. For compactness of the paper, only the modifications to apply the TO method in 3D are presented. For further details, readers are referred to Xia et al. (2020a) and the book by Bendsøe and Sigmund (2003). Compared to the TO procedure in 2D, the design structure (or full domain) is discretized by trilinear hexahedral structural finite elements for implementation convenience and numerical efficiency. The computational cost of 3D finite element method (FEM) analysis combined with iterative optimizations needs to be considered. In order to conduct the TO process within reasonable time, additional measures are taken. In assembling the global stiffness matrix K, because the components in the matrix are calculated independently and individually, the assembling process can be programmed in a parallel manner (Sharma & Martin, 2009), using multiple processing cores. For solving the governing equation (K(ρ)u(ρ) = f, where ρ, u, and f indicate element densities, nodal displacement, and nodal force, respectively), the preconditioned conjugate gradient method is adopted to iteratively solve this large, sparse linear system of equations. In this way, the required random access memory can be significantly reduced for 3D FEM analysis in comparison to matrix factorization that was applied in 2D. The efficiency of applying the iterative solvers in the TO process has been discussed by Amir and Sigmund (2011). Finally, the dimension of the global stiffness matrix K is reduced during the optimization process. The degrees of freedom (DOFs) that are solely linked to void elements (ρ i = 10 −4 ) are detected. A small number is used as the void density in TO to prevent singularity of the system matrix without affecting the TO result. The components K i associated with these detected DOFs are eliminated in K, which also speeds up the TO process. In order to verify the effectiveness of these measures, the computation time of analyzing 3D FEM models of various meshes is presented in Table 1. Without these measures, the computation time increases exponentially. Each measure improves the computation performance, and the most efficient procedure is obtained by adopting all of them. In the case of 30 × 30 × 30 elements, the computation time of adopting three measures is reduced to 3.5% of the original one, that is 27.81 over 784.56.

2.2
Step 2: TO extraction The obtained TO results (voxel meshes and densities) of Step 1 cannot be directly used for defining STM models. Postprocessing steps are usually required to convert TO results into required models for manufacturing purposes (Liu & Ma, 2016). For a review of 2D TO extraction techniques, the reader is referred to our previous paper (Xia F I G U R E 1 Flowchart of the 3D-OPT-STM generation method et al., 2020a). For 3D TO results, CAD-friendly models were extracted from TO results in Cuillière, François, and Nana (2018) and Yin, Xiao, and Cirak (2019). None of these methods is dedicated to, or suited for, obtaining STM models. In order to address this gap, in our previous study (Xia et al., 2020a), an automatic extraction method was developed to convert 2D TO results into truss-like models. In this section, the extraction method is extended and improved to generate truss-like models in the 3D setting. The TO extraction method includes three substeps (Step 2 in Figure 1). The obtained TO densities show intermediate values. By setting a threshold value (0.1), all densities larger than this value are set to 1 and the remainder is set to 0. In this way, the noise in the thinning process by the intermediate densities is avoided. The densities are transformed into binary data for the subsequent skeletonization step.

F I G U R E 2
The thinning process applied to a table-like model Based on the binary density and voxel mesh information, the raw TO results are simplified by a 3D digital thinning method. The thinning method reduces a large amount of TO data to a simple centerline pattern (skeleton) that preserves the topology. Based on the obtained skeleton, the subsequent truss-like model identification step can be implemented effectively. The pattern-based thinning method in our previous 2D work (Xia et al., 2020a) cannot be extended to 3D because the prescribed patterns are not able to describe all elimination conditions in 3D while preserving topology. In this paper, therefore the 3D thinning method developed by Lee, Kashyap, and Chu (1994) is adopted. The invariance of the Euler characteristic and the preservation of the connectivity are considered in 3D thinning. The Euler characteristic refers to the number of connected objects, holes, and cavities, which quantifies the topology of a given voxel model. The obtained skeleton must have the same Euler characteristic as the original voxel model to preserve the topology. The thinning method iteratively eliminates border voxels until obtaining a onevoxel width skeleton. A skeleton of a table-like model after the thinning process is shown in Figure 2. Note that the loaded and support points in our study are taken as unremovable voxels; they remain the same in generated STM models.
The truss-like model consists of a 3D network of nodes and connected bars. The extraction procedure in our previous 2D work is not suitable for 3D conditions. The previous extraction method may lead to overly complex 3D skeletons containing superfluous members. In this paper, a new procedure is proposed to extract 3D truss-like structures. In the procedure, the nodes of a 3D truss-like model are first identified, next their connections (i.e., the bars of the truss-like model) are determined, and finally short bars are removed to avoid difficulties in later STM shape optimization. Although these general steps are similar to the 2D version (Xia et al., 2020a), there are important differences in the details as described below.
The nodes of a 3D truss-like model are identified based on element-wise detection of skeleton voxels. Each 3D F I G U R E 3 Process of identifying nodes and their connections for the generation of a 3D truss-like model from a voxel-based skeleton skeleton voxel is surrounded by 26 neighboring voxels. If the number of solid neighboring voxels is greater than 2, this skeleton voxel is taken as a node (see in Figure 3a). Different from the procedure in our 2D work, clustered nodes are allowed in a small region. They will be removed in the following merging step while preserving the topology. After determining the nodes, the presence of direct connections via the skeleton indicates the bars of a truss-like model. In order to identify these connections, a 3D recursive detecting method is proposed because the complexity of 3D topologies makes our previous 2D connection identification process infeasible. The method node-wisely takes every identified node as the starting point for a connecting bar. By tracing the skeleton via the 26 neighboring voxels of each visited voxel, a path along the skeleton to another identified node can be found, which is taken as the bar end point (see in Figure 3b). Therefore, a bar of a trusslike model can be created by a straight line connecting the identified start and end nodes. The initial generated trusslike model may have various short bars that are impractical and insignificant for the STM model. Therefore, a merging method is proposed to eliminate these short bars. In the method, a minimum length is set to check all bars. If the bar length is smaller than the prescribed value, the short bar and its end nodes are replaced by a new node at its center location. Based on the proposed extraction method, a truss-like model of the table-like model ( Figure 2) is generated, as shown in Figure 4. The robustness of the proposed F I G U R E 4 The generated truss-like models generated by the proposed extraction process based on the table-like model (see in Figure 2) procedure is demonstrated by various case studies in the later sections.

2.3
Step 3: STM shape optimization The STM method is a truss-analogy based method: axial force equilibrium is necessary. The importance of the satisfaction of axial equilibrium has been demonstrated in Xia et al. (2020a). As the generated truss-like structures from the previous phase are usually statically and kinematically unstable structures, their equilibrium forces cannot be calculated through truss analysis (Xia et al., 2020b). In this paper, 3D beam elements with high slenderness (height/length = 10 −3 ) are used to calculate the equilibrium forces of the generated truss-like structures. Although equilibrium forces are obtained with beam analysis, the existence of non-marginal shear forces does not satisfy the requirement of the STM method. Thus the generated truss-like structure cannot be used as an STM without obtaining a pure axial force equilibrium. In order to obtain a suitable model, a 3D shape optimization method is introduced in this section.
In the shape optimization, the node positions of the generated truss-like structure are adjusted to optimize the structure. The mathematical formulation of this shape optimization is shown as where x is a vector containing the coordinates of the nodes of the generated truss-like structure, which are used as optimization variables. C, f, u, and K here represent the compliance, nodal forces, nodal displacements, and the stiffness matrix, respectively, based on the beam finite element model. There are two more considerations in the shape optimization: the axial equilibrium force state and the minimum bar length. They are considered through the constraints in the optimization process. Firstly, the suitable truss structure (STS) index is a dimensionless number that indicates the closeness of a beam structure to a pure axial force equilibrium state (Xia et al., 2020b). For a structure consisting of 3D beam elements, STS is given as where, is the element axial force, V e indicates the element shear forces in two orthogonal transverse directions, and n is the number of elements. STS = 1 indicates that a pure axial force equilibrium state is obtained for a specific beam structure. As shear forces are inevitable in the beam analysis, in this paper, a tolerance ɛ = 0.05 is used to relax the STS requirement in the optimization process. Second, in the inequality constraint ( ) ≥ min , min is the given minimum bar length to avoid short bars. L(x) measures the minimum length of all bars. Different from the previous 2D method, this constraint prevents the generation of short bars that would lead to convergence difficulties. Using a p-norm formulation, the measure of the minimum bar length is calculated as follows: where L e indicates the element length and n is the number of elements and the exponent is set to p = 10. The continuous minimum length constraint facilitates solving the STM shape optimization problem using gradient-based optimizers.
The shape optimization problem given in Equation (1) is solved using the gradient-based method of moving asymptotes (Svanberg, 1987). Central finite differences are used to calculate sensitivities. The perturbation is taken as 0.1% of the mesh size of the continuum elements. In order to ensure robust convergence of the nonconvex shape optimization problem, a move limit of 50% of the mesh size of the continuum elements is imposed on the optimization variables. Moreover, with the shape optimization, the structural compliance is further improved by correcting the differences caused by simplifying TO continua to discrete structures.

F I G U R E 5 Two simple cases and resulting structures after the 3D shape optimization
In order to briefly illustrate the effectiveness of the shape optimization, two simple cases are presented as shown in Figure 5. The first case (Figure 5a) is topologically a line with a helical shape under compression. The resulting STS index is 0.42. After the 3D shape optimization, as expected, a straight structure (Figure 5b) is obtained that carries the load by axial compression. The STS index of the optimized structure is 0.99 after 51 iterations. The second case (Figure 5c) is an irregular planar frame with two vertical forces and two hinges at two bottom points. It is a modification of a classical STM model (see Schlaich & Schafer, 1991) by shifting all members and adding redundant nodes. The STS index of this structure is 0.73. After 25 shape optimization iterations, the optimized model (Figure 5d) results in STS = 0.98 and it has a similar shape as the classical STM model. As the shape optimization involves a computationally inexpensive beam model, the whole optimization process takes less than a minute.

THREE-DIMENSIONAL OPTIMIZATION-BASED STRUT-AND-TIE MODELS FOR THREE TYPICAL D-REGIONS
In this section, the proposed generation method is applied to three typical D-regions: a four-pile cap, a corbel, and a box girder. Apart from the generated 3D-OPT-STM mod-F I G U R E 6 Geometry, load, and supports of the four-pile cap els, also three manually created STM models are considered. In all cases, the strength of the struts is verified and the strength of the ties is ensured by providing sufficient reinforcement in the direction of the ties. The capacities of the 3D nodal zones are considered as sufficient. In order to demonstrate the effectiveness of the generated 3D-OPT-STM models, the STM designs are simulated by NLFEA to obtain structural performances, in particular load-displacement curves and insight in structural failure modes. The structural efficiency is defined as the ratio between load capacity and steel usage. Finally, the steel usage, analyzed structural performances and structural efficiency between the generated 3D-OPT-STM models and manually created STM models are compared.

Four-pile cap case
Reinforced pile caps are widely used in bridges to transfer loads from piers to the pile foundation. They are common structures and previous studies have also proposed STM models for them. Here we first consider this case to verify that the proposed method generates a suitable and highquality STM model. Subsequently, in Section 4.1, with our systematic generation method, we explore STM models of a wide range of pile caps with different aspect ratios. The geometry and load conditions of the evaluated fourpile cap case are shown in Figure 6. In the figure, L, W, and H indicate the length, width, and height of the structure, respectively. In the considered model, the structure is subjected to a concentrated load of 700 kN and supported with four ball joints. Figure 7 illustrates the various steps of the generation method of obtaining the 3D-OPT-STM for this case. The FEM model of the design problem is presented in Figure 7a, with a mesh size of 12 mm and 62,500 solid elements. Concrete material properties are used: Young's modulus is 30 GPa and the Poisson ratio is 0.2. In the TO process, only linear finite element analysis is conducted. A volume constraint of 5% and a filter radius of 250 mm (2.5 × (mesh size)) are adopted. The binary TO result is  Figure 7b. The skeleton after the thinning process and the truss-like structure after the identification process are shown in Figure 7c,d, respectively. The obtained 3D-OPT-STM is shown in Figure 7e with an increase in STS index from 0.77 to 0.99, which is close to a pure axial force equilibrium state. In this case, using a PC with an Intel Core 4 3.5 GHz processor and 16 GB memory, the overall computational time is about 40 min. The 3D TO step uses about 95% of the total time, the extraction step is the second most expensive step and followed by the shape optimization, which takes less than 1 min, but offers crucial improvements that enable the use of the generated model in the STM method. Based on 3D-OPT-STM, the force diagram can be calculated as shown in Figure 7f. It can be used for the subsequent STM analysis.
In this paper, a manually created STM model is used for comparison, which is similar to the models in Dey and Karthik (2019), Souza et al. (2009), andYun et al. (2019), as shown in Figure 8. It has the same node positions as the generated 3D-OPT-STM, but a different topology, and results in an STS index of 0.98. Based on the obtained tensile forces, bar lengths, and steel yielding strength (580 MPa), the required steel according to the two STM models can be calculated.
In order to gain further insight into the structural performance of these two STM models, two NLFEA models F I G U R E 8 The manually created STM model are created and analyzed to obtain their structural performances. The NLFEA models of the two STM models, material properties of concrete and steel, and used solution strategies are described in the Appendix. Loaddisplacement curves of these two STM designs based on the NLFEA results are shown in Figure 9. In these, the displacement is measured vertically at the middle-bottom point of the cap. The 3D-OPT-STM model results in a peak load of 876.6 kN and the classic STM model leads to a peak load of 830.1 kN. Both peak loads are larger than the design load (700 kN). In the figure, Point A indicates the onset of cracking, and Points B and C indicate the steel yielding of the 3D-OPT-STM model and the classic STM model, respectively. It is observed that both structures fail due to steel rupturing.

F I G U R E 9 Load-displacement results of the 3D-OPT-STM design and classic STM design of the four-pile cap
Based on the generated STM models, the steel ratio of the generated 3D-OPT-STM and the manually created STM model are both equal to 0.39%. Although the topologies are different, they lead to the same steel ratio in this case. In Blévot and Frémy (1967) and Clarke (1973), the similar performance of these two designs was observed in their experiments. By comparing the ratio of peak load and steel usage, 2.081 and 1.971 N/mm 3 for the 3D-OPT-STM and classic model, respectively, more efficient steel usage is observed in the 3D-OPT-STM design. In a later section (Section 4.1), further investigations on applying the generation method for the four-pile cap with various geometrical parameters are presented.

Corbel case
Reinforced concrete corbels are common D-regions used in building structures to support floors or precast beams.
Here we first verify the proposed method to generate a suitable STM model for this case. In addition to the plane load conditions, corbels may be loaded in out-of-plane direction, which will be studied in Section 4.2 as well. For this paper, the FEM model and structural conditions of the evaluated corbel are shown in Figure 10. The corbel is fixed at the top and bottom ends and a concentrated load of 300 kN is applied. The FEM model consists of 14,337 solid elements with a mesh size of 33 mm. Similar to the fourpile cap case, concrete material properties are used.
In the TO process, the corbel is optimized with a 10% volume constraint and a 49.5 mm filter radius (1.5 × (mesh size)). After 107 TO iterations, the TO result of this corbel case is shown in Figure 11a. The resulting 3D-OPT-STM F I G U R E 1 0 The FEM model of the corbel case (mm) F I G U R E 1 1 Generation results of the corbel case. Red and blue bars indicate tensile and compressive members, respectively. The bar width indicates the force magnitude is obtained after 26 shape optimization iterations and is shown in Figure 11b. The STS index has increased from 0.84 to 0.97. Based on the calculated tensile forces and steel yielding strength (580 MPa), the required steel ratio of this corbel is 0.18%. The overall computational time of this case is about 10 min. Similar to the four-pile cap case, the TO step accounts for 95% of the computational time.
In Figure 12a, an STM model based on a 2D model by Schlaich et al. (1987) is manually created for this corbel case. This manually created STM model is formed by a stable truss and results in an STS index of 1. The calculated equilibrium forces are presented in Figure 12b, and this design leads to a required steel ratio of 0.22%. Compared to the manually created STM model, the 3D-OPT-STM design has a steel ratio of 0.18% and thus results in a reduction of 0.04% (percentage points of steel ratio) of steel usage and offers a more economical design. This corresponds to a relative reduction of 18%. Corresponding NLFEA models of these two STM models are presented in the Appendix. The resulting load-displacement curves are shown in Figure 13. The displacement is measured at the middle-bottom point of the outer end of the corbel. Based on the NLFEA results, the peak loads of the 3D-OPT-STM model and the classic STM model are 336.1 and 378.0 kN, respectively. The ultimate capacity of both designs is larger than the design load (300 kN). The concrete starts cracking at Point A shown in Figure 13. The steel yields at Point B for the classic STM model, and at Point C for the 3D-OPT-STM model. Both designs exhibit steel rupturing in the final stage. Although the classic STM design results in a larger peak load, by comparing the ratio of peak load and steel usage, 0.342 and 0.319 N/mm 3 for the 3D-OPT-STM and classic model, respectively, a more efficient steel usage is observed in the 3D-OPT-STM design.

Box girder case
Box girders are often used in bridge structures. Compared to common I-shaped or H-shaped beams, box girders provide larger load resistance capability, especially in resisting torsion. There is no standard STM model available for this fully 3D case. First, the effectiveness of the proposed method is verified through the generated STM model. Subsequently, the influence of complex load combinations (in Section 4.2) and load discretization (in Section 4.3) for the STM models are explored for the box girder based on the proposed method. The FEM model of an evaluated box girder is shown in Figure 14. The box girder is subjected to nine concentrated vertical forces of 300 kN, which represent a discretization of a uniformly distributed load on the top flange. The structure is constrained in axial direction at the end surface and is supported additionally by two ball joints. The FEM model has 118,400 solid elements with a mesh size of 50 mm and has the same concrete material properties as the four-pile cap case.
In the TO process, the box girder is optimized with a 10% volume constraint and a 125 mm filter radius (2.5 × (mesh size)). After 134 TO iterations, the TO result of this corbel case is shown in Figure 15a. The resulting 3D-OPT-STM is obtained after 57 iterations of the shape optimization and is shown in Figure 15b. The optimization procedure leads F I G U R E 1 5 Generation results of the box girder case. Red and blue bars indicate tensile and compressive members, respectively. The bar width indicates the force magnitude F I G U R E 1 6 The manually created STM of the box girder. Red and blue bars indicate tensile and compressive members, respectively. The bar width indicates the force magnitude to an increase in the STS index from 0.77 to 0.98. Based on the calculated tensile forces and steel yielding strength (580 MPa), the required steel ratio of this box girder design is 0.175%. The overall computational time of this case is about 250 min. In this case, the TO step accounts for about 95% of the computational time.
As the manually created STM model for this box girder, we use the spatial truss structure as shown in Figure 16a. The manually created STM model results in an STS index of 0.97. The calculated equilibrium forces are shown in Figure 16b. This STM design leads to a steel ratio of 0.245%. Compared to the manually created STM design, the 3D-OPT-STM design thus resulted in a reduction of 0.070% of steel usage (percentage points of steel ratio). This corresponds to a relative reduction of 29%.
Based on the NLFEA models (as shown in the Appendix) for the two STM models, load-displacement curves are presented in Figure 17, where the displacement is measured at the middle-bottom point at the free edge. In the curves, Point A indicates the initial cracking of concrete and Points B and C indicate the steel yielding points of the 3D-OPT-STM design and the manually created STM design, respectively. Both models reach the design load level of 2,700 kN. After obtaining the design load, the models indicate a rapid failure due to steel rupturing. The man-F I G U R E 1 7 Load-displacement results of the 3D-OPT-STM design and classic STM design of the box girder ually created STM model results in a slightly larger peak load, however substantially more steel is required. The ratios of peak load to steel usage of the 3D-OPT-STM and manual STM models are 0.114 and 0.0775 N/mm 3 , respectively. Also in this case, the 3D-OPT-STM design produces a more economical design than the manually created STM design. With these findings, the generation method can be considered capable of generating high-quality STM designs, in a practical timeframe.

THREE IMPORTANT FACTORS IN GENERATING 3D STM MODELS
The generated STM models will vary due to changes in input parameters and problem setting, such as the geometry of a specific case, the load conditions and the load discretization. Based on the proposed generation method, these factors can be investigated in a convenient and systematic manner. Here this capability is demonstrated for three different 3D cases. First, various 3D-OPT-STM models are generated by varying geometrical parameters of the four-pile cap. The steel usage of the generated 3D-OPT-STM models is compared with that of classical STM models. Second, the load conditions influence TO results and consequently affect the generated 3D-OPT-STM models. Generating 3D-OPT-STM models considering alternative load conditions is investigated. Third, load discretization is a prerequisite of the STM and affects the generated STM models. Therefore, the influence of the load discretization for the generated 3D-OPT-STM models is discussed.
F I G U R E 1 8 Generated 3D-OPT-STM models of various geometries. The gray voxels in each model indicate TO results. L, W, and H indicate the length, width, and height of the four-pile cap, respectively. The numbers in the gray boxes indicate the steel ratio, in percentage. Right of those, the green/black/red numbers indicate the difference (also in percentage) compared to STM models with a classical topology given in Figure 8

Parametric investigation for the four-pile cap
Four-pile caps with various geometries are used in practice. In this section, 25 3D-OPT-STM models are generated by adjusting the width and height of the four-pile cap, represented by (W, H) in Figure 6a. In the various models the length (L = 600 mm) is fixed. The width and height are modified accordingly using five variations, as H/L = [2/6, 3/6, 4/6, 6/6, 12/6] and W/L = [2/6, 3/6, 4/6, 5/6, 6/6], respectively. By varying these parameters, the span-depth ratio and plane shape of the cap are changed. The parametric settings of the evaluated case in Section 3.1 are H/L = 3/6 and W/L = 6/6.
Following the same generation process with the fourpile cap case (Section 3.1), the obtained 3D-OPT-STM models are shown in Figure 18. The required steel ratios for all 3D-OPT-STM designs are given. Accordingly, classical STM models are evaluated, which have the same topology as the model shown in Figure 8. Com-paring the required steel ratio of 3D-OPT-STM designs with classical STM designs, the differences are shown in Figure 18. The green and red numbers indicate the reductions and increases of steel ratios by applying the automatically generated 3D-OPT-STM designs, respectively.
Based on the calculated steel ratios, we observe that by increasing the H/L ratio, which indicates a reduction of the span-depth ratio of the structure, the steel ratios reduce. By increasing H the internal level arm becomes larger, making the longitudinal reinforcement more efficient. For very large values of H, the results show a direct load transfer between the loaded point and supports. By reducing the W/L ratio, the load transfer is approaching a one-way system (i.e., beam action). Although the steel ratio increased, due to the reduction of the structural volume the total steel usage is reduced. Compared to classical STM models, the 3D-OPT-STM models reduce the steel usage in the situation of the larger span-depth ratios. In cases where the 3D-OPT-STM models have higher steel ratios, the F I G U R E 1 9 The corbel problem under complex load combination. The gray voxels indicate the TO results. Red and blue bars indicate tensile and compressive members, respectively. The bar width indicates the force magnitude differences are very small. In most cases, the 3D-OPT-STM models have a similar performance in terms of steel usage as the classical STM models.
Based on the obtained results, the best scenario in reducing steel usage is obtained in the model of (W/L, H/L) = (4/6, 2/6). Compared to the classical STM model, the 3D-OPT-STM design results in a reduction of the steel ratio of 0.349% (percentage points of steel ratio). Especially for low and wide pile caps, using the STM models generated by our method instead of the classical one leads to significant steel savings. Note also that the proposed generation method enables performing a parametric investigation in an automatic and convenient manner.

The influence of complex load combinations
In practice, D-regions are subjected to various loads. The suitable STM model depends on the considered load combinations. Especially in the 3D setting classical 2D STM models are definitely not applicable for cases with loads in out-of-plane directions. Clearly, in such cases considering the full 3D situation is necessary.
In order to investigate the influence of the complex load combinations, two cases are studied. The first case is based on the corbel problem (Section 3.2), and a horizontal load is considered, as shown in Figure 19a. This load may come from the supported structures, such as bridge girders and crane beams. They transfer the horizontal load to the corbel. Using the same settings of the generation method as in the previous case (Section 3.2), the 3D-OPT-STM model is obtained as shown in Figure 19b. Compared to the 3D-OPT-STM model with one single load (Figure 11b), the model considering combined loads leads to an obviously different load transfer system. Although the obtained 3D-OPT-STM model has an irregular shape, it results in an STS of 0.97. The twisted shape of the obtained model is beneficial in resisting torsion caused by the horizontal force. The steel ratio based on this 3D-OPT-STM design is 0.239% which is smaller than the single-load design (0.28%). Apparently, with the right design less steel is required in the case with multiple loads, for this problem. However, more loads are transferred to the supports through compressive forces in struts. In this case, the corbel would fail by strut crushing under the combined load situation.
The second case is similar to the box girder problem (Section 3.3), as shown in Figure 20a. Here, two compressive forces are considered, which represent applied prestressed tendons. The prestressed tendons are commonly used in box girders to improve structural service performance, such as reducing cracks and deformation. Similarly, based on the same generation settings with the singleload problem, the generated 3D-OPT-STM model for the combined loads is shown in Figure 20b.
The generated 3D-OPT-STM has a 0.98 STS index. Compared to the single-load 3D-OPT-STM model (Figure 15b), the load transfer system is different. Here the generated model has fewer tensile members. The loads are mainly transferred by compression through struts to the supports. The steel ratio based on the generated model is 0.048%, which is much lower than that of the single-load model (0.175%). By applying the post-tensioned tendons, the required steel in the box girder is significantly reduced. However, if we account the steel usage of the prestressed tendon with the same yielding strength (580 MPa), the additional steel ratio of the tendon is 0.175% and the total steel ratio is 0.223%. This can be further improved by changing the position of prestressed tendons.
Based on the previous two cases, the effectiveness of the proposed generation method for complex combined loads has been demonstrated. The generated 3D-OPT-STM models under combined loads are different from the single-load models. In both cases the amount of reinforcement reduces, however this is not a general rule. The locations of steel reinforcement are clearly influenced by the presence of additional loads. Moreover, in these two cases, more compressive forces are obtained compared to the single-load situation. Therefore, the strut strength checking in the STM analysis becomes even more important.
In practice, many loads and load combinations must be dealt with. In these cases, various 3D-OPT-STM models need to be generated. The proposed generation method facilitates the process in finding and exploring suitable STM models under various load situations, in limited time (up to a few hours for detailed models) and in an automated manner. Note that in the studied cases, all loads are assumed to occur simultaneously. For less predictable loads, such as earthquake or wind loads, a multi-load case formulation must be considered, which is an extension left for future research.

The influence of load discretization
Discrete loads are required in the STM. Different load discretizations affect the generated STM models and thus their performances. Especially in the 3D setting, the spatial distribution of discrete loads may have a significant influence on the generated STM models. In order to investigate the influence of load discretization choices, three 3D-OPT-STM models of the box girder case (Section 3.2) are generated. In Section 3.2, the uniform vertical load is uniformly discretized into 3 × 3 concentrated loads of 300 kN. Now two more load discretization schemes are applied, 2 × 2 loads of 675 kN and 4 × 4 loads of 168.75 kN. The obtained TO results and 3D-OPT-STM models are shown in Figure 21.
The generated two 3D-OPT-STM models result in STS indices of 0.97. Based on the obtained tensile forces, the required steel ratio of three load schemes (2 × 2, 3 × 3, and 4 × 4) are 0.165%, 0.175%, and 0.160%, respectively. In F I G U R E 2 1 Generated 3D-OPT-STM models of two additional load discretization schemes. The gray voxels indicate the TO results. Red and blue bars indicate tensile and compressive members, respectively. The bar width indicates the force magnitude the 2 × 2 scheme, the loads are aligned with the two webs of the box and the two supports. This results in a nearly 2D STM model. Similarly, the 2D-like performance is also observed in the model of the 4 × 4 scheme. The 3D-OPT-STM model of 3 × 3 scheme results in a slightly larger steel usage than the other two cases: additional steel in the top flange is required to transfer the load from the central loading points to both webs. The main tensile forces of three cases occur in the bottom flange, 1,605, 1,593, and 1,513 kN, respectively, and their differences are small. However, the location and forces of the main compressive members are changed in these load conditions.
The proposed generation method enables fast and systematic exploration of the influence of various load discretization schemes. A suitable load discretization scheme is required for a safe and robust STM design. In this box girder case, a more robust result would be obtained through the 3D-OPT-STM model with 3 × 3 discrete loads. It presents a more refined steel distribution on the top flange to resist the possible bending behavior.

CONCLUSION
In this paper, we propose a systematic method to automatically generate STM models in 3D. The proposed method integrates a TO process, the topology extraction method and the shape optimization procedure, to automatically generate suitable 3D-OPT-STM models for the design of 3D D-regions. In the TO, three measures (parallel matrix assembly, iterative solution techniques, and removal of redundant DOFs) are combined to improve the computation efficiency. A new robust procedure is proposed to extract truss-like structures from 3D TO results. Subsequently, gradientbased 3D shape optimization is conducted to effectively obtain suitable 3D-OPT-STM models based on the extracted truss-like structures. Depending on the complexity of the model, within one or several hours an OPT-STM model is produced. Based on the proposed generation method, three 3D-OPT-STM models are generated for three typical 3D D-regions (four-pile cap, corbel, and box girder) to demonstrate the effectiveness of the proposed method. Moreover, based on the generation method, three important factors in STM are discussed, including the parametric investigation for the four-pile cap, the influence of different load combinations, and the load discretization. Based on the present investigation, the conclusion is summarized as follows: 1. The proposed method is effective in generating suitable 3D-OPT-STM models for 3D D-regions in an automatic manner. The generated and used STM models are simulated by NLFEA, and they are evaluated based on steel usage and the simulated structural performance. Compared to manually created STM models, the 3D-OPT-STM models were found to produce economically superior designs, especially for complex problems. It is reiterated that practical aspects, such as steel detailing and constructability were outside the scope of our study. 2. Parametric investigations can be conveniently conducted by the proposed method. It enables engineers to evaluate a larger range of design options, even at an early design stage. The change of the load transfer mechanism is observed through the generated STM models for different geometries. 3. The proposed method enables the systematic generation of valid 3D-OPT-STM models under different load conditions. Based on the obtained 3D-OPT-STM models, various load transfer systems are presented which can help the designer to conduct a more accurate STM analysis. In this study, all loads are assumed to occur simultaneously. In order to generate more robust and safer STM models, a multi-load case formulation is of interest for future research. 4. Load discretization schemes are important in 3D STM analysis and affect the generated 3D-OPT-STM models. A more robust design can be obtained by choosing a suitable load discretization. The suitability of various load discretization schemes can be readily and systematically studied by the proposed method.

APPENDIX A NONLINEAR FINITE ELEMENT ANALYSES OF STRUT-AND-TIE MODELS
In NLFEA, nonlinear phenomena, such as the cracking and crushing of concrete and the yielding and rupturing of steel, are considered. The simulation results include the loading capacities and failure modes, which are used to evaluate the STM models.
In the NLFEA models, we assume that reinforcement bars are located corresponding to the ties in the STM model and their cross-sections are based on full utilization of the yielding stress at the design load. A similar solution strategy of NLFEA is used in our previous work (Xia et al., 2020a) for 2D STM models. The embedded truss elements are used for rebar, assuming perfect bonding with concrete. The Hordijk model and the parabolic model are used to present the tensile and compressive behavior of concrete, respectively, and the rotating smeared crack model is adopted to indicate the crack initiation (Belletti, Damoni, & Hendriks, 2011). In this paper, 20-node quadratic solid elements are used. The concrete has a mean compressive strength of 30 MPa and its derived properties are given in Table A1. The steel has a yield strength of 580 MPa with- The NLFEA model of the 3D-OPT-STM design for the four-pile cap is shown in Figure A1. A nodal displacement of 0.1 mm for a single-load step is applied vertically on the top steel plate. Four steel plates are used to support the cap. For the steel plates, a linear material model is adopted. Similarly, the analysis model of the manually created STM model was created by modifying the rebar while the concrete remained unchanged.
The NLFEA model of the 3D-OPT-STM design for the corbel is shown in Figure A2. A nodal displacement of 0.1 mm for a single-load step is applied vertically on the steel plate. The corbel is fixed at the top and bottom surfaces. The model for the classic STM design is created by modifying the rebar.
The NLFEA model of the 3D-OPT-STM design for the box girder is shown in Figure A3. In this case, a vertical equally distributed load of 0.3 MPa is applied to the nine surfaces (200 mm × 200 mm) on the top. A similar NLFEA model is created for the manually created STM model by modifying the rebar.