Innovative Approach for Enhancing GLULAM Performance with Reinforcing Steel Bars: A BESO-based Study

Abstract Glued-Laminated Timber (GLULAM) is a widely-used building material, popular for its strength, durability, and sustainability. It is created by bonding together layers of wood, making it a common choice for civil structures. In this study, Bi-directional Evolutionary Structural Optimization (BESO) is proposed to improve the performance of GLULAM structures. By positioning steel bars within the GLULAM structure, the objective is to increase the structure's stiffness and enhance its structural integrity. To achieve this, the study introduces the concept of a sub-design domain and utilizes optimization theory to determine the optimal placement of the steel bars. The finite element problem is solved using ANSYS software, while the topological optimization problem is solved using MATLAB software. The use of sub-design domains and optimization theory enables the optimal placement of the reinforcements to be determined. The results of this study demonstrate the potential of this approach for enhancing the structural integrity and stiffness of GLULAM structures under static loads. Graphical Abstract


INTRODUCTION
Engineered wood products, such as Glued Laminated Timber (GLULAM) and Cross Laminated Timber (CLT), are gaining popularity as substitutes for conventional construction materials like concrete and steel.GLULAM is created by stacking and bonding wood layers in the same direction as the wood grain.At the same time, CLT involves rotating the grain direction relative to neighboring layers by approximately 90 degrees.These engineered wood products are widely used in the construction industry in countries such as the USA, Canada, Europe, and Asia, due to their increased use in civil construction (Kremer & Symmons, 2015;Lall et al., 2019).Although the cost per volume of GLULAM and CLT is higher compared to steel and concrete, it is a low-density material that offers cost savings in transportation and foundation.Additionally, engineered wood structures are typically transported in semi-ready forms, which helps keep the workspace clean and saves time in the assembly process (Gustavsson et al., 2006).
The increasing demand for carbon reduction (D'Amico et al., 2021;Hart et al., 2021) and the advantages of engineered wood products have driven studies to enhance their structural capacity (Franzoni et al., 2016(Franzoni et al., , 2017(Franzoni et al., , 2018;;Peixoto et al., 2021).The integration of engineered wood products with other materials is being explored as a strategy for decarbonizing new buildings (D'Amico et al., 2021).This can be achieved by combining wood with materials such as concrete (Sebastian et al., 2017), aluminum (Chybiński & Polus, 2019), bamboo (Sun et al., 2020), and steel (Chiniforush et al., 2018;Hassanieh et al., 2016aHassanieh et al., , 2016bHassanieh et al., , 2017;;Peixoto et al., 2021;Soriano et al., 2016).Soriano et al. (2016) increased the overall mechanical capacity of GLULAM by incorporating 10 mm diameter steel bars through grooves fabricated symmetrically along the neutral axis and the entire length of the GLULAM.Substituting traditional civil construction materials with wooden structures like GLULAM and CLT has the potential to reduce CO2 emissions by up to 31% (Oliver et al., 2014).By 2050, the use of low-carbon substitutes for specific applications could prevent the emission of 50 Mt CO2 in certain regions (D'Amico et al., 2021;Hart et al., 2021).
Exploring slight modifications in the geometry and size of traditionally manufactured timber structures is a continuous field of study.For instance, small gaps in CLT structures can lead to a loss of mechanical strength.However, the use of intentional gaps filled with other materials in engineered wood products has the potential to drive innovation in design and improve acoustic and thermal properties (Franzoni et al., 2016(Franzoni et al., , 2017(Franzoni et al., , 2018)).
Using topology optimization tools can lead to improved design features that meet both the model's constraints and user requirements during the design process.This optimization method is used to determine the optimal material distribution within a defined region, and there are several approaches, including the solid isotropic material with penalization (SIMP) method (Bendsoe & Sigmund, 2004;Zhou & Rozvany, 1991), the evolutionary structural optimization (ESO) method, the bi-directional evolutionary structural optimization (BESO) method (Huang & Xie, 2007;Xie & Steven, 1993), and the level set method (Sethian & Wiegmann, 2000).
Several areas of academia use topological optimization as a subject of study or as a tool, including minimizing the frequency responses of multiscale systems (Vicente et al., 2016a) composed of macro and micro phases (Vicente, Zuo, et al., 2016b;Wang et al., 2022), multi-material topology optimization methods considering isotropic and anisotropic materials and their combination (Bohrer & Kim, 2021), and topology optimization considering fluid-structure iteration (Picelli et al., 2015;Vicente et al., 2015).
In engineered wood products, the optimization method was applied in a truss topology optimization framework involving two materials, wood and steel (Ching & Carstensen, 2022).The Global Warning Potential (GWP) was proposed, indicating steel and wood's global warming capacity.When allied with compliance, it is possible to find truss structures with stiffness compromise and low GWP involving wood and steel (Ching & Carstensen, 2022).Mayencourt and Mueller (2019) suggested a novel topological optimization method for determining the ideal relative density of the intermediate layers in one-way slab CLT structures.They found that the optimal relative density was 0.41 for a 5-layer CLT, leading to the use of alternating solid wood and void spaces to achieve this density.The implementation of this method resulted in a significant decrease in material consumption, reducing costs by 18% with only a minimal loss in capacity.
In addition to the potential for material and cost reduction, CLT with an optimized core may have advantages in acoustic insulation when used as flooring.Huang et al. (2021) experimentally evaluated the core optimization methodology of CLT (Mayencourt & Mueller, 2019) when subjected to human-generated sound pressure.The authors found that CLT with an optimized core did not worsen the sound insulation capacity and, in some cases, had higher capacity than CLT without voids.Perković et al. (2021) evaluated GLULAM structures with hollow cores of elliptical and circular shapes and found that these structures can experience a loss of up to 40% of their structural capacity.de Vito et al. (2023) applied the BESO method to GLULAM and CLT structures, aiming to minimize displacement under static force conditions while maintaining a constrained volume in 3D models.Their proposed model considers the orthotropic properties of wood and the stacking of layers.However, the authors found that not implementing periodicity constraints in topological optimization may negatively impact the fabrication of optimized structures.
Latin American Journal of Solids and Structures, 2023, 20(6), e503 3/13 To the best of author's knowledge, the literature has yet to investigate the topological optimization of timber structures with reinforcement such as steel rods (Mayencourt & Mueller, 2019, 2020;Peixoto et al., 2021;Soriano et al., 2016).In this sense, this work aims to develop a methodology based on the BESO method to find the best positioning of reinforcement for GLULAM structures subjected to static loads and volume constraints.For this purpose, this work is divided into sections, as follows: Section 2, entitled 'Methods,' presents the developed method, including the definition of the topology optimization problem, material property, material interpolation, sensitivity analysis, filter, optimization criteria, and strategy for using the sub-design domain.Section 3 entitled Numerical Implementation,' presents 2D and 3D numerical examples illustrating the methodology presented in the previous section.Finally, in Section 4, the conclusion of the study is presented, along with future work and observations derived from the completion of this research.

METHODS
The method used to apply topological optimization with multiple materials in this paper is based on the consolidated BESO methods (He et al., 2016;Huang & Xie, 2007, 2008).For maximizing the stiffness of the structure composed of two materials, the mean compliance should be minimized according to the problem stated in Eq. 1: where u is the displacement matrix, and  is the stiffness matrix; F and u are the values of external forces and displacements applied in the structure domain.Ku = F stands for the equilibrium equation of the system.N stands for the number of elements in design-domain,  1 * is the fraction of the prescribed volume for the steel beam and  2 * is the fraction of volume for the wood.∑    =1 stands for total design domain volume.  stands for element density of th element for th material ( = 1 steel and j = 2 wood), Given that if  = 1, then   = 1, and if  = 2,   =   .

Material and Material Interpolation
The cases evaluated in this paper involve two solid phases: steel and wood.Young's modulus of wood is considered orthotropic according to values taken Forest Service & Products Laboratory (2010) for the Douglas fir species in the longitudinal, radial, and tangential axes.The torsional strength and Poisson's ratio are also determined using the same method.Figure 1 represents the main axis of the wood in relation to the fiber direction, Longitudinal, Radial and Tangential.Steel is treated as isotropic.The parameters used in this work are presented in Table 1 for reference.
When dealing with topology optimization problems that involve multiple phases, determining the material distribution between these phases is often accomplished using an interpolation equation that describes the relationship between them.
Table 1 Assumed mechanical properties.

Material Parameters
For two solid phases such as steel and wood, the material interpolation can be expressed using the penalty factor p in the following equation:

Sensitivity Analysis
The sensitivity analysis determines the importance of each element in the design-domain.The sensitivity number of the compliance due to an element change in the design domain is obtained by deriving the objective function for the design variable

∂𝐶𝐶 ∂𝑥𝑥 𝑖𝑖𝑖𝑖
. This term can be obtained by differentiating the equilibrium equation: Since changing an element from solid to void does not change the applied force, the term ()   = 0. Appling the chain rule in Eq. 3: Deriving the equilibrium equation gives that: Applying Eq. 4 in Eq. 5, it follows that: Whereas, from equilibrium equation    − =   the sensitivity number is written as: Latin American Journal of Solids and Structures, 2023, 20(6), e503 5/13 Based on the presence of two solid phases and the material interpolation described in Eq. 2, the sensitivity number is derived by quantifying the change in mean compliance or total strain energy resulting from the modification of an element in a structure (X.Huang & Xie, 2008) .Thus, the expression for the sensitivity number is as follows: where in the context of two solid phases and the material interpolation presented in Eq. 2, the stiffness matrix for the ℎ element is calculated separately for isotropic steel �   � and orthotropic wood �   � characteristics.

Sensitivity Filtering, Stabilization of Evolutionary Process and Convergence Criterion
The sensitivity number must be filtered to avoid checkerboard-like problems.There are several effective methods presented in the literature (Chen et al., 2019;Zuo & Xie, 2015).In this paper, the filter applied is related to a filter size   determined by the user and defined by Eq. 9 and Eq.10: where α  is the filtered sensitivity number, the parameters  and   are dependent on the filter radius   .Where  is the number of elements within the Ψ sub-domain defined by a radius   , according to Figure 2; α  are the sensitivity numbers referring to the elements belonging to the Ψ sub-domain;   is the weight parameter defined in Eq. 8. where   is the centroid distance from element  to element .In addition to the checkerboard-like problems, the convergence can be improved by averaging the sensitivities values from the previous iteration (X.Huang & Xie, 2010; X. Huang & Xie, 2007), as described in Eq. 11: where α  � are the averaged sensitivities and  is the value of the current iteration.When the volume of the current iteration (  ) reaches the final volume (  * ) the optimization continues, without changing the volume for the next iteration, until the convergence criterion, described in Eq. 12, is reached.
where  is the convergence error. states the number of iterations of stable compliance set to 4.

Sub-design Domain Definition, Element Change Criterion and Update Scheme
The position of the steel bars within the wood is determined by discretizing the sub-design domain according to the specific requirements and availability of each problem.In this study, for the 2D example, the wood beam was divided into design domains of reinforcement beam sizes Ls and Hs, as shown in Figure 3.To optimize the positioning of the steel bars in the sub-design domains, it is necessary to calculate and filter the sensitivity number of each element in the structure according to Eqs. 8, 9, 10, and 11 at each iteration.Next, the elements belonging to each sub-design domain are identified.This allows for the adjustment of the sensitivity number of each element belonging to the sub-design domain by the average sensitivity number of that domain.In this way, the importance of each sub-design domain can be measured by the average of its constituent elements.The quantity of bars to be used in each iteration must then be defined.Assuming that the design-domain starts completely with steel ( = 1 and   = 1), in each iteration, a value ER is interpreted as the number of bars changed per iteration, which is used to calculate the number of bars in the next iteration.The sub-design domains with lower sensitivity numbers have their elements changed to wood  = 2 and   =   , obeying the number of bars altered per iteration.The simulation continues until the number of steel bars/number of sub-design domains is reached and until Eq. 12 is satisfied.Similarly, in 3D cases, the sub-design domain is defined by characterizing the steel bar size to be used, but with an additional dimension, Ws.In this study, steel bars with diameter d were used and positioned as shown in Figure 4.As with the 2D case, it is necessary to separate the elements belonging to the sub-design domain and rank the sensitivity number.However, in this case, only the elements belonging to the d diameter bar of the sub-design domain with the highest sensitivity number should be replaced with steel.This paper utilizes two software programs, MATLAB and ANSYS, to apply the proposed method.ANSYS is responsible for solving the finite element problem, while MATLAB is responsible for updating the parameters and solving the topology optimization problem.The sequence and indication of the software responsible for each step are shown in Figure 5.In the computational modeling of the wood composites, the orthotropic characteristic of the wood is considered.In addition, the glue between the layers of glued laminated timber structures was not considered, thereby being a rigid joint between the layers.The entire strategy involving finite element modeling in ANSYS software and calculation sequence in MATLAB is available at https://github.com/arturvito/Reinforced-GLULAM.
The sequence for applying the proposed method is adapted from (X.Huang & Xie, 2007).Before adding, removing, or altering a beam in the wooden structure, it is necessary to evaluate the fraction of steel reinforcement beams � 1 +1 � that will be required for the next iteration (k+1), where k represents the current iteration.Since the fraction specified by the user for reinforcement beams ( 1 * ) can be greater or smaller than the initial design, the volume of the next iteration is calculated in each iteration as follows: where ER is referred to as the evolutionary volume ratio, as presented by (X.Huang & Xie, 2007).However, in this work, it directly represents the quantity of reinforcements altered per iteration.Additionally, since the void phase is not considered, the value of  2 +1 is automatically calculated.Once the specified fraction of reinforcement beam volume is reached, the number of beams in the structure is not changed.With the calculated volume, it is necessary to calculate the sensitivity number for all elements of the structure using Eq.8.Unlike the approach proposed by (X.Huang & Xie, 2007), this work adds a step to the strategy, which is to identify all elements belonging to the sub-design domain and replace the sensitivity value with the average sensitivity value of all elements in that sub-design domain.The remaining steps follow the proposal by (X.Huang & Xie, 2007).Therefore, this work follows the following evolutionary iteration procedure: 1. Discretize the computational model through finite element meshing and boundary conditions.Define the elements belonging to the initial structure as j=1 steel and j=2 wood.
2. Solve the computational model using the finite element method and export the values for calculating the sensitivity number of each element.
3. Perform filtering and averaging procedures on the sensitivity numbers.
4. Identify the elements belonging to each sub-design domain and replace the sensitivity number values of each element with the average sensitivity number of the elements belonging to the sub-design domain.
5. Determine the quantity of reinforcement beams required for the next iteration.
6. Update the structure by adding, removing, or altering reinforcement beams.
7. Repeat steps 2-7 until the desired number of beams is achieved and the convergence criterion is satisfied.

NUMERICAL IMPLEMENTATIONS
This section presents two numerical examples of the method proposed in this paper.The first example is a 2D model where the initial topology consists of a GLULAM beam reinforced with steel in all possible positions, and steel structures are removed iteratively according to the stipulated methodology and parameters.The second example is a 3D application with a sub design domain composed of wood and a cylindrical bar in the center.These examples are intended to demonstrate the effectiveness and versatility of the proposed method for both 2D and 3D models.

Numerical Example 1 -2D GLULAM Beam Reinforced with Steel
In this example, only half of the model is simulated to decrease computational cost.The dimensions of the GLULAM are L = 2400 mm and H = 120 mm, supported at the end under a platform and subjected to a concentrated load in the center of the beam, as shown in Figure 6.The steel beam in this half of the model has dimensions of   = 150 mm by   = 6 mm.Initially, all 40 sub-design domains are filled with steel bars, and at each iteration, two bars are removed until a total of 8 bar segments remain.The other parameters for the BESO method are   = 65 mm, τ = 0.01, and p = 3.As the number of reinforcement bars decreases from 14, there is a significant increase in compliance.However, when the structure has more than 14 bars, compliance decreases at a gradual gradient.Moreover, for comparison purposes, a simulation was carried out using completely wooden conditions and results in a compliance value of 1.43 × 10 −3 Nm, which is approximately 68% lower than with the placement of 14 steel bar segments, indicating that the structure is 68% less resistant.Interestingly, this value increases to approximately 80% when 40 steel bar segments are used.
Figure 8 presents the distribution of steel bars obtained from four iterations of the simulation.The first iteration of the simulation represents the initial guess of the structure, which consisted of 40 steel bars.As the simulation progresses, steel bars are incrementally removed from the structure based on the filtered sensitivity number.In the eighth iteration, 16 steel bars were removed near the support and force application regions.In the 12th and 25th iterations, the bars were removed from the center, leaving steel bars at the ends of the structure near the fixed point and force application location.

Numerical Example 2 -3D GLULAM Beam Reinforced with steel
In this example, the method proposed in the paper is applied to simulate only half of the GLULAM model in order to reduce computational costs.The dimensions of the GLULAM were based on and adapted from the study by (Soriano et al., 2016).Each sub-design domain has a length of   = 200.0mm, height of   = 30.0mm, and width of   = 30.0mm.To create the 120mm x 120mm x 1400mm GLULAM half-beam finite element method model, four sub-design domains are required in both width and height, and seven sub-design domains in length, making a total of 112 possible locations for the bars.A 6mm diameter cylinder is positioned in the center of each sub-design domain, where it is defined as either wood or steel, while the rest of the sub-design domain is wood.The positioning of the cylinders and the finite element mesh used, half of the beam, and the applied load is shown in Figure 9.
Initially, 112 bars of steel are present in the structure, the steel bars are removed two-by-two until 28 bars remains.The other parameters for the BESO method were   = 30 , τ = 0.001, and p = 3.
Figure 10 displays the progress of the objective function for each iteration and steel bar variation.It is worth noting that, due to the limited number of steel bars, the simulation may reach a point where no further improvement in the objective function can be achieved.At the beginning of the optimization process, the removal of bars has a smaller influence compared to the end of the simulation.This observation suggests that an excessive number of steel bars in the GLULAM beam does not necessarily lead to a proportionate increase in the structure's rigidity.However, after 34 bars, there is an abrupt change in compliance, indicating a strong correlation between this number of bars and the structural strength.Figure 11 shows four possible configurations of the arrangement of 28 steel bar segments in the structure proposed in numerical example 3. The first configuration is the result of the applied method, which corresponds to a compliance of 322.35 Nm with the steel bars arranged in the central part of the beam and close to the fixed end.The next option presented is condition A, in which the 28 steel bar segments are arranged at the ends of the beam, resulting in a compliance of 346.56 Nm, approximately 6% higher than the optimized result.In condition B, the beams are arranged in the center of the width and at the ends of the height, resulting in a compliance of 346.16 Nm.Conditions A and B present similar compliance values, as there are beams positioned in the optimal regions and, therefore, are close to the best case for this condition.In condition C, the 28 steel bar segments are located in the center of both width and height, resulting in a compliance of 669.53 Nm, approximately 51% higher than the optimal condition.In this case, none of the steel bar segments are located in the optimal position.

CONCLUSION AND FUTURE WORKS
This study has demonstrated the effectiveness of utilizing the BESO method to optimize the positioning of steel reinforcing bars in GLULAM structures.By introducing the concept of sub-design domains and leveraging optimization theory, the proposed approach has enabled the determination of the optimal areas/volumes within the GLULAM structures to place the steel bars.The results of the simulations have shown that adding reinforcement bars to the wooden structure can significantly enhance its stiffness by up to 68%, exhibiting nonlinear behavior.However, the study also indicates that there is a threshold point at which the addition of bars reaches a plateau, highlighting the importance of carefully selecting the number and position of reinforcing bars.
The proposed method provides a valuable contribution to the field of engineering and construction by introducing a novel approach to optimizing GLULAM structures.The findings have demonstrated the potential for significantly enhancing the structural integrity and stiffness of GLULAM structures under static loads, which can lead to improved safety, durability, and cost-effectiveness.Furthermore, the use of GLULAM in construction reduces carbon dioxide emissions during transportation, making it a sustainable choice for building materials.
The results presented have important implications for the design and manufacturing of composite structures and could lead to more efficient and sustainable engineering solutions.Further research into other optimization techniques and their effectiveness for different types of composite materials and structural designs could provide valuable insights for the development of more effective optimization algorithms.Despite the simplifying assumptions made in the study, the authors acknowledge the practical aspects of GLULAM fabrication and recognize the need for further investigation into these optimized structures from a manufacturing perspective in future research.The proposed methodology can be extended to other types of structures and materials, making it a valuable tool for engineers and researchers in the field of structural optimization.

Figure 2
Figure 2 Schematic representation of filter sub-domain and radius   .

Figure 3
Figure 3 An example of the design domain and sub-design domain for 2D case.

Figure 4
Figure4An example of the design domain and sub-design domain for 3D case.

Figure 5
Figure 5 Sequence adopted for solving the topological optimization problem.

Figure 6
Figure 6 Design domain and sub-design domain of example 1 with indication of concentrated load application and boundary conditions.

Figure 7
Figure7provides valuable insight into the optimization process, illustrating how the objective function changes as the simulation progresses.Specifically, the graph shows the relationship between the objective function and the number of reinforcement bars used in the structure.

Figure 7
Figure 7 Optimization histories of the objective function and reinforcement bar quantity for example 1.

Figure 8
Figure 8 Four intermediate topologies with indication of the steel bars position from numerical example 1.

Figure 9
Figure 9 Design domain and sub-design domain of example 3D GLULAM Beam Reinforced with steel with indication of distributed load application region and boundary conditions.

Figure 10
Figure 10 Optimization histories of the objective function and reinforcement bar quantity for example 2.

Figure 11
Figure 11 four possible configurations for the positioning of steel bar segments in the GLULAM structure.From left to right, the first configuration is the result of topological optimization.The next configuration, condition A, has the steel bar segments located at the ends of the GLULAM beam.In condition B, the bars are positioned at the top and bottom ends of the GLULAM beam.Finally, in condition C, the bars are located in the center of the GLULAM beam.