Numerical Evaluation of Lateral Torsional Buckling of PFRP Channel Beams under Pure Bending

Abstract

The use of pultruded fiber reinforced polymers (PFRPs) in strengthening and sustainable design of bridges and other structures exposed to corrosion and resistance reduction factors is growing rapidly. However, a comprehensive understanding of the structural behavior of these materials under various loading conditions is crucial to unlock their full potential and promote their wider use in diverse structural and industrial applications. Pultrusion profiles can be also used as beams in bridges. One important aspect of the structural behavior of PFRPs is their buckling behavior, particularly in thin-walled open cross sections. Lateral torsional buckling is a probable instability mode for beams with thin-walled open cross sections that are not laterally restrained along their span. Therefore, research on the buckling behavior of PFRP members is essential. In this study, the analytical responses of channel-shaped PFRP beams in bridges under pure bending are calculated using an equation in the Eurocode 3 regulation. The buckling behavior of these beams is then investigated through numerical modeling using the finite element package Abaqus. A total of 75 specimens of PFRP channel profiles with different thicknesses in various spans and lateral restraint conditions are studied for their lateral-torsional buckling behavior. This study uniquely explores the behavior of PFRP beams with lateral restraints, a novel aspect in the field of lateral-torsional buckling research of PFRP beams. The results show that the analytical equation used for these beams needs to be modified to more accurately estimate the buckling loads of FRP beams under the conditions studied in this paper.

1. Introduction

Sustainability is a critical consideration in the design and maintenance of bridges, particularly in relation to corrosion and the overall lifespan of the structure. Bridges play a vital role in transportation networks and are essential for economic development and societal well-being. However, the deterioration of bridge infrastructure due to corrosion and other factors poses significant challenges in terms of safety, maintenance costs, and environmental impact.

Corrosion, caused by factors such as exposure to chloride ions and aggressive environments, is a major concern for bridge structures. It leads to the deterioration of materials, reduces the structural integrity of bridges, and necessitates costly repairs and maintenance to extend their service life [1]. Therefore, sustainable design and maintenance strategies that address corrosion prevention and mitigation are crucial for ensuring the long-term viability of bridge infrastructure.

In recent years, the integration of sustainability into the design and upkeep of bridges has gained significant attention. Various research efforts have been directed at making bridge structures more sustainable by focusing on key aspects like energy efficiency, safety standards, corrosion onset, and impact on the environment [2]. A vital element in this pursuit is the role of structural health monitoring (SHM) and the reliability of bridges, which are both fundamental to their sustainable development. Research such as [3] has made strides in utilizing sophisticated Bayesian networks for the analysis of structural system reliability. This work improves the precision and effectiveness of reliability evaluations, an essential factor in the sustainable oversight of bridge infrastructure. Furthermore, the study in [4] brings a novel perspective by applying sparse Bayesian learning to SHM data, leading to a more accurate data modeling approach that is critical for dependable and long-lasting SHM systems. Adding to this, the research presented in [5] highlights the advantages of using probabilistic and data-centric methods for detecting structural damages. This approach is key to enhancing the predictive maintenance and extending the lifespan of bridge structures. Moreover, the work in [6] introduces an innovative method for modeling and predicting strain measurements in bridges, especially during severe weather conditions like typhoons. This contribution is notably crucial for improving the resilience aspect of sustainable bridge engineering. Together, these studies mark a significant advancement in the realms of structural health monitoring and reliability analysis, guiding the sustainable evolution of bridge infrastructure. Additionally, the use of multi-objective decision-making methods has proven effective in evaluating and choosing the most sustainable design options [7,8]. Such techniques allow for a comprehensive consideration of various factors and uncertainties, enabling decisions that are well-informed and aligned with sustainability goals.

Furthermore, the use of innovative materials and technologies has emerged as a promising avenue for sustainable bridge design. For example, the application of fiber reinforced polymers (FRPs) in bridge strengthening and construction offers corrosion resistance and enhanced durability [9]. Among the various manufacturing techniques, pultrusion stands out as a process capable of producing continuous FRP profiles. The utilization of these profiles, particularly Pultruded Fiber Reinforced Polymer (PFRP) segments in construction, even as the main structural elements of buildings [10], has been steadily increasing due to their advantageous characteristics, including lightweight design, high resistance to fatigue, freedom from corrosion, and favorable bending strength [11,12]. PFRPs have shown potential for improving the structural behavior of bridges and reducing the need for frequent maintenance. Understanding the behavior of these materials under various loading conditions, including buckling behavior, is crucial for their effective utilization in bridge design [9].

Composite beams can undergo various types of buckling, including local buckling, global buckling involving lateral, torsional, or flexural modes, as well as instabilities resulting from the interaction between local and global buckling. Extensive research has been conducted to investigate the local buckling behavior of FRP beams, encompassing different cross-sectional profiles [13,14,15,16]. These investigations have employed diverse methodologies, encompassing experimental investigations [17,18,19,20] as well as analytical and computational approaches [21,22,23]. Similarly, the global buckling phenomenon has been the subject of multiple studies, exploring buckling related to flexure, torsion, and lateral force individually or in combination, across a range of section types such as I, angle, rectangular, channel, and strip sections.

One of the critical stability concerns for beams with thin-walled open cross-sections arises when they lack sufficient lateral restraint along their span, leading to potential lateral-torsional buckling. On the other hand, the design considerations for a PFRP member are predominantly influenced by elastic deflections and/or elastic buckling instabilities rather than strength, due to its relatively low stiffness-to-strength ratio [24,25]. Consequently, it is crucial to focus on studying the buckling behavior of PFRP materials in order to gain valuable insights. This paper is dedicated to investigating the structural behavior of PFRP profiles, which are characterized by channel-shaped cross sections and are fabricated using FRP composite materials via the pultrusion method. The primary objectives of this research encompass assessing the lateral-torsional buckling strength of PFRP profiles and evaluating their deformations when subjected to pure bending loading. By addressing these key aspects, valuable insights can be gained to enhance the understanding and application of PFRP profiles in structural engineering.

This study presents a unique investigation into the behavior of PFRP beams with lateral restraints, a topic that has not been extensively explored in existing research. While previous studies have delved into various aspects of lateral-torsional buckling in PFRP beams, including the effects of unrestrained lengths and the buckling behavior under different loading conditions, the specific focus on the influence and control of LTB through lateral restraints represents a novel contribution of our work. This approach provides new insights into enhancing the stability of PFRP beams, making a significant contribution to the field of sustainable bridge and structural design.

2. Numerical Modeling

The finite element method (FEM) has been employed to simulate the response of these beams using the Abaqus finite element package. The composite beams were modeled based on the information and characteristics provided by a manufacturer of these profiles, wherein E-type glass fibers and polyester resin were used as the reinforcing material and polymer, respectively [26]. In this study, macroscopic material properties, drawn from the accessible data, are utilized for modeling the samples. In earlier research, like Nguyen’s research [27,28] and the experimental study done by Zeinali et al. [29], this method has received validation, yielding results that align closely with those derived from laboratory experiments.

2.1. Materials and Geometrical Properties

Based on experimental investigations, it has been determined that the longitudinal strains observed during the lateral-torsional buckling (LTB) failure of a PFRP beam remain in the elastic range [13,30]. Consequently, the material modeling approach employed is based on linear elasticity. It is essential to incorporate orthotropic properties when modeling PFRP materials due to their varying mechanical characteristics across three mutually perpendicular directions. These directions, denoted as 1, 2, and 3, correspond to the fiber’s longitudinal direction, the transverse direction, and the height of the sample, respectively. As a result, the mechanical properties of the materials are specified in the longitudinal, transverse, and altitudinal directions to account for these perpendicular orientations.

Table 1. Material properties.

${\mathit{E}}_{\mathit{L}}$ (MPa)	${\mathit{E}}_{\mathit{T}}$ (MPa)	${\mathit{\mu }}_{\mathit{L}\mathit{T}}$	${\mathit{\mu }}_{\mathit{T}\mathit{T}}$	${\mathit{G}}_{\mathit{L}\mathit{T}}$ (MPa)
23,000	8500	0.23	0.09	3000

provides an overview of the physical properties of PFRP materials. In numerical modeling, nine engineering constants are utilized to define orthotropic PFRP materials. These constants include $E_1,E_2=E_3,G_{23},G_{12}=G_{13}$, and ${\mu }_{23},{\mu }_{12}={\mu }_{13}$, where $E_1,E_2$, and $E_3$ represent the three moduli present in orthotropic materials. The PFRP materials in the beams possess a layered arrangement, resulting in identical elastic moduli in the transverse and altitudinal directions. Specifically, $E_1$ corresponds to the longitudinal modulus of elasticity $\left(E_L\right)$, while $E_2$ and $E_3$ represent the transverse modulus of elasticity $\left(E_T\right)$. The longitudinal and transverse elastic moduli of the beam are determined by the combined elastic moduli of the fibers and resin. Additionally, $G_{23},G_{12}$, and $G_{13}$ denote the three shear moduli in orthotropic materials, which are equivalent to the in-plane shear modulus $\left(G_{LT}\right)$. Parameters ${\mu }_{23},{\mu }_{12}$, and ${\mu }_{13}$ represent three Poisson’s ratios in orthotropic materials, which are equal to the major Poisson’s ratio (${\mu }_{LT}$). Based on specifications provided by the manufacturing company and experiments conducted on these materials [31], the difference in shear modulus across different directions is less than 10%. Hence, it is assumed that the shear modulus values are uniform across all three directions. These specifications serve as input values for defining the materials within the modeling framework.

For the purpose of modeling in this study, a comprehensive selection of 75 specimens was made, each featuring three distinct channel sections characterized by dimensions of 120 × 50 × 10 mm${}^3$, 120 × 50 × 8 mm${}^3$, and 120 × 50 × 6 mm${}^3$ (height × width × thickness). These specimens were carefully chosen to represent a range of spans, employing simply supported boundary conditions. The considered spans for the beams ranged from 1.8 to 4.2 m. To achieve varying slenderness ratios, specific lateral restrained distances were implemented on the upper flange, which experiences compressive forces. Ultimately, the specimens were simulated under pure bending loading to accurately analyze their structural behavior.

Table 2. Samples with a thickness of 6 mm.

	$\mathit{L}/{\mathit{L}}_{\mathit{b}}$	1	1.2	1.6	1.8	2
$\mathit{L}/\mathit{D}$
15		S1 = 1800	S2 = 1500	S3 = 1080	S4 = 1000	S5 = 900
20		S6 = 2400	S7 = 2000	S8 = 1440	S9 = 1300	S10 = 1200
25		S11 = 3000	S12 = 2500	S13 = 1800	S14 = 1670	S15 = 1500
30		S16 = 3600	S17 = 3000	S18 = 2160	S19 = 2000	S20 = 1800
35		S21 = 4200	S22 = 3500	S23 = 2520	S24 = 2300	S25 = 2100

,

Table 3. Samples with a thickness of 8 mm.

	$\mathit{L}/{\mathit{L}}_{\mathit{b}}$	1	1.2	1.6	1.8	2
$\mathit{L}/\mathit{D}$
15		S26 = 1800	S27 = 1500	S28 = 1080	S29 = 1000	S30 = 900
20		S31 = 2400	S32 = 2000	S33 = 1440	S34 = 1300	S35 = 1200
25		S36 = 3000	S37 = 2500	S38 = 1800	S39 = 1670	S40 = 1500
30		S41 = 3600	S42 = 3000	S43 = 2160	S44 = 2000	S45 = 1800
35		S46 = 4200	S47 = 3500	S48 = 2520	S49 = 2300	S50 = 2100

and

Table 4. Samples with a thickness of 10 mm.

	$\mathit{L}/{\mathit{L}}_{\mathit{b}}$	1	1.2	1.6	1.8	2
$\mathit{L}/\mathit{D}$
15		S51 = 1800	S52 = 1500	S53 = 1080	S54 = 1000	S55 = 900
20		S56 = 2400	S57 = 2000	S58 = 1440	S59 = 1300	S60 = 1200
25		S61 = 3000	S62 = 2500	S63 = 1800	S64 = 1670	S65 = 1500
30		S66 = 3600	S67 = 3000	S68 = 2160	S69 = 2000	S70 = 1800
35		S71 = 4200	S72 = 3500	S73 = 2520	S74 = 2300	S75 = 2100

classify the beams according to various parameters, including their length (L), restraint distance ($L_b$), length-to-height ratio ($L/D$), and beam length-to-restrained length ratio ($L/L_b$). The values in the tables represent the restrained lengths of the beams, measured in millimeters. For instance, all samples in the first row have a total length of 1800 mm, but the restrained length of S2 is 1500 mm. Similarly, all samples in the second row have the total length of 2400 mm.

2.2. Loading and Boundary Conditions

In order to establish simply supported conditions, it is imperative to define the coordinate directions for the specimens and subsequently establish the appropriate boundary and support conditions with respect to these directions. For the numerical models, the longitudinal, transverse, and height directions are aligned with the X, Y, and Z axes, respectively. Consequently, the displacements in the X, Y, and Z directions within the models are denoted as $U_X,U_Y$, and $U_Z$, respectively. Similarly, the rotations around these axes are represented as $U{R}_X,U{R}_Y$, and $U{R}_Z$, respectively.

To induce pure bending loading, the beams are subjected to loading at the shear center of the section in such a manner that the upper flange experiences compression, while the lower flange experiences tension. In Figure 1, the definition of shear center using a reference point and rigid body in numerical modeling is shown. Adhering to the predefined directions, the boundary conditions are set as $U_X=U_Y=U_Z=0$ at one support and $U_X=U_Y=0$ at the other support (as shown in Figure 2). The location of the support coincides with a transverse line on the lower flange. Additionally, lateral restraints are applied exclusively to the upper flange ($U_Y=0$), along with the specific distances. These restraints restrict movement solely in the lateral direction. Furthermore, lateral restraints are implemented at the support locations to prevent the beam from overturning on the upper flange.

2.3. Element Types and Mesh Sizes

The beams in this study were modeled using shell elements, which are well-suited for representing low-thickness planar structures. Shell elements disregard the analysis through the thickness of the structure. Among the suitable element types for investigating the buckling behavior of PFRP beams, S4R, S4R5, and S8R elements were considered [27]. S4R and S4R5 elements were found to be more suitable for thin shells, while the S8R element was better suited for thick shells. In this study, the S8R shell element type is employed for modeling purposes (Figure 3). The S8R element is an 8-node, quadrilateral shell element that represents a stress/displacement doubly curved thick shell. It incorporates reduced integration and utilizes a large-strain formulation [32].

The modeling approach utilized integrated meshing, with the dimensions of the mesh influencing the sensitivity analysis. Figure 4 illustrates the mesh sensitivity analysis curve for the load-lateral displacement diagram of one of the samples. The finer mesh has dimensions of 5 and 10 mm, while the coarser mesh has dimensions of 25 and 50 mm. Comparing these figures reveals that increasing the number of meshes results in coinciding load-lateral displacement diagrams. A mesh size of 12.5 mm is chosen, with four elements present in each cross-sectional flange (Figure 5), representing the deformation of a sinusoidal wave, if it occurs. Finer meshing below 12.5 mm significantly increased the problem-solving time but had no noticeable effect on the results. Similar observations were made for other samples.

2.4. Analysis Type

The modeling process involves two steps. In the first step, eigenvalue analysis (linear buckling analysis) is employed to determine the critical loads for elastic buckling in a pure bending state, specifically for the first mode. During eigenvalue analysis, the beam is analyzed in the absence of any imperfections. Moving on to the second step, the beam is modeled with a slight curvature, serving as an initial imperfection located at the midpoint of the beam span. An imperfection of 1/1000th of the sample length is introduced. Nonlinear buckling analysis, specifically the Riks method [33], is then employed to determine the critical loads for elastic buckling in a pure bending state.

3. Closed Form Equation

As previously mentioned, pure bending is often the most critical loading mode for beams. Under this loading condition, the likelihood of LTB is higher compared to other loading states, and the LTB strength is lower. According to Eurocode 3, the torsional buckling moment ($M_{cr}$) for a steel beam with symmetry about the major axis is determined using Equation (1) [34]. For PFRP materials, the critical buckling moment can be calculated by modifying Equation (1) and replacing the physical properties of steel with the corresponding properties of PFRP materials [35,36,37,38].

$$M_{cr}=C_1\frac{{\pi }^{2}E{I}_z}{{\left(kl\right)}^2}\left\{{\left[{\left(\frac{k}{k_w}\right)}^2\frac{I_w}{I_z}+\frac{{\left(kL\right)}^{2}G{I}_t}{{\pi }^{2}E{I}_z}+{\left[C_2{z}_g\right]}^2\right]}^{0.5}-C_2{z}_g\right\}$$

(1)

Equation (1) includes various parameters: E represents the modulus of elasticity, G denotes the shear modulus, k refers to the effective length factor against lateral bending, and $k_w$ represents the equivalent factor for end warping. In the case of pure bending, both k and $k_w$ are assumed to be equal to 1. The span length of the simply supported beam is represented by L. Additionally, Equation (1) incorporates $C_1$, an equivalent uniform moment factor that accounts for the shape of the bending moment distribution. In pure bending, $C_1$ is considered to be equal to 1. $C_2$ is a factor that considers the vertical load height relative to the shear center, with $Z_g$ representing the height of the load from the shear center. When the load is positioned above the shear center (towards the top flange), $Z_g$ is positive, and when it is below (towards the bottom flange), $Z_g$ is negative. In this study, since we are considering the pure bending state, both $C_2$ and $Z_g$ are assumed to be zero.

The second moments of area for flexure about the minor axis ($I_Z$), warping rigidity ($I_W$), and torsional rigidity ($I_t$) are denoted in Equation (1). For channel sections used in this study, the values corresponding to $I_Z$, $I_t$, and $I_W$ can be obtained from

Table 5. Geometric properties of the sections.

Section Dimensions (mm)	${\mathit{I}}_{\mathit{Z}}$ (mm${}^4$)	${\mathit{I}}_{\mathit{t}}$ (mm${}^4$)	${\mathit{I}}_{\mathit{W}}$ (mm${}^6$)
$120\times 50\times 6$	$2.78\times {10}^5$	$1.4\times {10}^4$	$6.2\times {10}^8$
$120\times 50\times 8$	$3.5\times {10}^5$	$3.5\times {10}^4$	$7.6\times {10}^8$
$120\times 50\times 10$	$4.2\times {10}^5$	$6.7\times {10}^4$	$8.6\times {10}^8$

and used in Equation (1).

4. Results and Discussion

Through a comprehensive analysis of the samples, the buckling load values obtained from Equation (1), as well as the results from numerical buckling analysis, are carefully examined. Additionally, load-displacement diagrams associated with the midpoint of the free span of beams are meticulously compared through nonlinear analysis. This comparative assessment aims to gain insights into the agreement and discrepancies between the different approaches, with a particular focus on the impact of lateral restraints. This research, among the first in its field, provide a comprehensive understanding of the buckling behavior exhibited by the beams under investigation.

4.1. Linear Buckling Analysis and Theoretical Examination

In the context of numerical linear buckling analysis, the software output typically includes the critical load value and the corresponding buckling mode shape. Other output information may not be provided in this step. Figure 6 demonstrates the buckling of the compression flange in a beam without any restraint, and the buckling shape is observed in another sample with restraint.

Based on the analysis of the moment-length-to-height ratio diagram derived from Equation (1) shown in Figure 7, it can be observed that the buckling moment exhibits a nonlinear decrease as the ratio of beam length to section height increases. Additionally, the effect of section thickness on the increase of the buckling moment diminishes as the three graphs converge towards the end. For instance, when the ratio of length to height of the section increases by 100% (from 15 to 30), the buckling moment decreases by approximately 60% for a section thickness of 6 mm, around 57% for a thickness of 8 mm, and approximately 54% for a thickness of 10 mm.

When considering buckling analysis, it becomes evident that the results obtained from linear analysis closely align with the values derived from the equation when the beam’s length equals that of the unrestrained span; in other words, the beam in the compression flange is without restraint. It shows that the linear analysis has reliable results when it is used for beams with span length but without lateral restraint. Another point regarding the results of linear analysis is that, as the lateral restraint approaches the roller support, the buckling load resulting from the linear analysis of the beam increases. The reason for this issue can be the creation of some kind of rigidity in the beam and, as a result, the buckling load increases. In general, the linear analysis of buckling does not provide us with the correct results of the buckling load. This type of analysis in this research is done only to introduce an imperfection in the structure to perform a more detailed Riks analysis. Of course, it should be mentioned that these statements are for the difference between the closed form equation and the numerical buckling analysis of the beams; however, another analysis has been done for these beams, which has investigated the nonlinear behavior of the beams, as discussed in Section 4.2.

4.2. The Results of Nonlinear Riks Analysis

In nonlinear analysis, the load is incrementally applied in a static manner until the sample either reaches buckling or yields. Once buckling initiates, the load remains relatively constant, while the displacement of the model significantly increases or the sample becomes completely unstable, causing the load to decrease. This type of analysis, unlike linear buckling analysis, is particularly valuable in cases where the presence of restraints makes it challenging to accurately model the buckling behavior. Nonlinear analysis provides more reliable and interpretable results, enhancing the accuracy of the predicted buckling response. Figure 8 shows the final deformation of samples S1 and S2, which are with and without lateral restraint, respectively.

In this paper, the load-lateral displacement diagrams are defined as the main output of the software. The purpose of including these diagrams is to show the difference in the buckling capacity of the samples and the lateral displacement in the mid-span of the compression flange of the channel-shaped beams. The vertical axis of the diagrams is related to the load, and its unit is Newton. The horizontal axis shows the lateral displacement in the middle of the free span of the compression flange, with a unit of millimeters. The sizes and specifications of the samples are those of Table 2, Table 3, Table 4 and Table 5.

4.2.1. Comparison of Samples with Different Length-to-Restrained-Length Ratios

Figure 9 is the load-lateral displacement diagram for the first five samples. The samples are 1.8 m long, 6 mm cross-sectional thickness, cross-section length-to-depth ratio of 15, with different length-to-restrained-length ratios (from 1 to 2), which have become unstable due to lateral-torsional buckling. As depicted in Figure 9, the sample labeled as S5, with the lateral restraint positioned at the midpoint of the beam, exhibits the highest buckling load among the analyzed samples. This finding suggests that the location of the lateral restraint significantly influences the buckling behavior. Additionally, it was observed that, as the critical length decreases, the buckling load increases. This equation highlights the importance of considering the critical length in predicting and understanding the buckling response of the beam.

As can be seen, the beam deviates from the initial path by bearing the load at one point and changes its direction to the positive side of the graph. This point of separation of the diagram from the initial path can be introduced as the beginning of lateral buckling-torsion in the load-lateral displacement diagram. In fact, the load in which the displacement change is significant compared to the previous state is considered the buckling load. According to the examination of the samples, it can be said that, with the increase in the ratio of the length to the length of the restrained span of the sample, the amount of buckling load has increased nonlinearly.

According to Figure 9, the presence of restraint has caused the buckling load to be at least 1.2 times and at most 3.3 times the buckling load of the corresponding sample without restraint. The minimum increase is related to the case where the restraint is closest to the support, and its maximum is the case where the restraint is placed in the middle of the span.

Another point is that only the critical length of the samples is considered in the regulations; however, by comparing the nonlinear analysis of samples with the same critical length and different span length, it was observed that the closer the beam span is to the critical length, the higher the buckling load. For example, an examination of

Table 6. Results of samples with thickness 6 mm.

Sample Name	Buckling Load from Linear Analysis (N)	Buckling Load from Nonlinear Analysis (N)	Buckling Load from Equation (1) (N)	Difference between Linear Response and Equation (1) (%)	Difference between Nonlinear Response and Equation (1) (%)
S1	1781.9	1703	2634.3	47.8	55
S2	2371.9	2650	3488.3	47	32
S3	2476.4	3765	6019	143	60
S4	2445.5	3942	6884.9	181	75
S5	2459.8	4150	8304	237.6	100
S6	1284.5	1183	1748.8	36	48
S7	1628.6	1770	2256.6	38.6	27
S8	1793.8	2530	3722.9	107	47
S9	1745.5	2596	4397.8	152	69
S10	1768.3	2711	5026.9	184	85
S11	1010.5	970	1306.6	29	35
S12	1240.2	1379	1655	33	20
S13	1457.4	1850	2634	81	42
S14	1393	1943	2950.7	112	52
S15	1424	2034	3488	145	71
S16	836.13	788	1044.5	25	32
S17	1004.9	1000	1306.6	30	31
S18	1267.1	1475	2022	59.6	37
S19	1185.6	1600	2256.6	90	41
S20	1224.8	2400	2634	115	9.8
S21	799.62	655	871.6	9	33
S22	847.53	880	1080	27	23
S23	1148.9	1267	1637.8	42.5	29
S24	1053.7	1920	1853	75.9	−3
S25	1094.5	1976	2104	92	6

,

Table 7. Results of samples with thickness 8 mm.

Sample Name	Buckling Load from Linear Analysis (N)	Buckling Load from Nonlinear Analysis (N)	Buckling Load from Equation (1) (N)	Difference between Linear Response and Equation (1) (%)	Difference between Nonlinear Response and Equation (1) (%)
S26	2949.1	2637	2949	26	49
S27	3679.2	3979	3679	32	27
S28	4049.5	5647	4049	77	47
S29	3943.8	5926	3944	107	58
S30	3993.6	6134	3994	138	82
S31	2174.3	2123	2174	21	28
S32	2618.9	2652	2619	27	29
S33	3148.5	3869	3148	44	39
S34	2993.1	3986	2993	80	56
S35	3068.3	4243	3068	93	66
S36	1729.4	1537	1729	17	36
S37	2044.3	2151	2044	23	21
S38	2704.3	2988	2704	22	32
S39	2491	2989	2491	54	46
S40	2593.1	4700	2593	66	7.7
S41	1438.1	1358	1438	16	25
S42	1683.9	1768	1684	21	18
S43	2451.7	2494	2452	18	25
S44	2161.9	2513	2162	44	36
S45	2288.5	2576	2288	52	53
S46	1232.1	1175	1232	15	22
S47	1435.4	2075	1435	20	−15
S48	2157	2792	2157	15	−7.8
S49	1913.2	3073	1913	40	−6.4
S50	2040.2	3157	2040	48	2

and

Table 8. Results of samples with thickness 10 mm.

Sample Name	Buckling Load from Linear Analysis (N)	Buckling Load from Nonlinear Analysis (N)	Buckling Load from Equation (1) (N)	Difference between Linear Response and Equation (1) (%)	Difference between Nonlinear Response and Equation (1) (%)
S51	4385.6	3950	5545	26	40
S52	5274.4	5753	6983.8	32	21
S53	6201.2	8261	10,974	76.9	33
S54	5937.8	8290	12,292	107	48
S55	6065.5	8478	14,425.9	137.8	70
S56	3269.7	3209	3943	20.6	23
S57	3847.6	3984	4880	26.8	22
S58	5097.5	5633	7367	44	31
S59	4697.8	5783	8451	79.8	46
S60	4889.9	6014	9441	93	57
S61	2612.7	1986	3071.8	17	55
S62	3046.6	2987	3764	23	26
S63	4549.9	4486	5545	21.8	24
S64	3955.3	4565	6087	53.9	33
S65	4203.2	6553	6983.8	66	6
S66	2177.4	2162	2521	15.8	17
S67	2529.7	2651	3071.8	21	16
S68	3786.5	3843	4455	17.6	16
S69	3389	3855	4880	44	26
S70	3634	5510	5545	52	0.6
S71	1867	1764	2141	14.6	21
S72	2165.8	3000	2598.7	19.9	−13
S73	3227.4	4095	3730	15.5	−8.9
S74	2959.4	4553	4141	39.9	−9
S75	3118.7	4690	4605	47.6	−1.8

reveals that sample S2 has a span length of 1.8 m, sample S15 has a span length of 3 m, and both have a restrained length of 1.5 m. Other geometric characteristics of both samples are also similar. This occurs when the buckling load in sample S2 is about 30% higher than the other sample, while the buckling load resulting from Equation (1) is the same for both samples.

By carefully comparing the buckling load of other samples, such as S27 and S40 or samples S52 and S65, an increase of 18% and 14%, respectively, in the buckling load has been observed. It can be said that, with the increase in thickness, the effect of the difference between the length of the span and the restrained length on the buckling load decreases.

4.2.2. Beams with Different Length-to-Section-Depth Ratios

In this section, "load-lateral displacement" graphs related to a group of five samples are drawn. Samples in this group of diagrams have the same section thickness and same length-to-restrained-length ratio, with different lengths having different length-to-section-depth ratios. All of the specimens are unstable due to lateral-torsional buckling.

Figure 10 is related to beams with a section thickness of 6 mm, with length-to-depth ratios of 15 to 35 and a length-to-restrained-length ratio of 1. In fact, these samples do not have lateral restraints along their span. In these samples, by increasing the ratio of length to depth of a section from 15 to 35, the buckling load is reduced by 160%, which corresponds to sample S21 compared to sample S1.

It is unnecessary to include repeated diagrams here, but by observing similar examples, it can be said that beams with a section thickness of 8 mm and 10 mm with a length-to-restrained-length ratio of 1 exhibit a reduction in buckling load by a maximum of 125% and 124%, respectively.

According to the outputs of this section, it can be said that the buckling load decreases with the increase in the length-to-depth ratio of the cross-section of the samples. It can also be said that the effect of the increase in the ratio of the length to the depth of the cross section of the samples on the buckling load decreases with the increase in thickness.

4.2.3. Output Diagrams of Beams with Different Section Thicknesses

In this section, a group of samples that have the same length, length-to-depth ratio, length-to-restrained-length ratio, and different section thicknesses is compared.

Figure 11 shows the diagrams of “load-lateral displacement” for beams with a length of 1.8 m, a ratio of the length to the restrained length of 1, and a ratio of length to depth of section equal to 15. These samples do not have lateral restraints along their opening. The thickness of sample 26 is 33% more than the first sample, while its buckling load is about 55% more than the first sample, and the 25% increase in thickness has increased the buckling load by 50%.

By comparing the same groups, we can see that, with the increase in the thickness of the pultruded FRP beams, the buckling load has increased nonlinearly. Another point is that, by increasing the ratio of the length to the length of the restrained span, the effect of increasing the thickness on the buckling load decreases.

4.2.4. Comparison of the Results from Different Analyses and Conclusions

According to Table 6, Table 7 and Table 8, the difference between the results of the linear analysis and the equation is 9% in the best case and 237.6% in the worst case, and the difference between the results of the nonlinear analysis and the equation is 0.6% in the best case and 100% in the worst case.

Figure 12 is drawn for a better comparison of the nonlinear analysis response and the response obtained from Equation (1). This diagram is related to a group in which S1, S6, S11, S16, and S21 are selected as examples. By looking carefully at the diagram, it can be seen that, with the increase in the ratio of the length to the height of the beam section, the buckling load resulting from Equation (1) and the nonlinear response are close to each other. Considering the constancy of the cross-section height in all samples, it can be said that the longer the beam length, the smaller the effect of the orthotropic behavior of the beam on the buckling load and the closer to the isotropic behavior of the beam. These results apply to other groups as well and, due to non-repetition, other graphs have not been included.

Figure 13 illustrates the effect of restrained distance on the difference between the results of nonlinear analysis and Equation (1). The isotropic equation provides a closer answer in case of beams with a greater restraining distance. It can be said that the isotropic equation for orthotropic elastic beams has a more accurate response and the orthotropic behavior of beams has a greater effect on the buckling load at less restraining distances. This diagram is drawn for beams with lateral restraint with a ratio of length to height of section 15, thickness of section 6 mm, and different restraining distances.

The significant difference between the results of Equation (1) and the nonlinear analysis indicates that the equation used is not accurate enough for these materials. There is a need to modify the equation that can more accurately estimate the buckling load of FRP beams investigated in this research. Considering the characteristics of the materials of these beams, which are considered orthotropic, and the equation used for isotropic materials, it can be said that a large percentage of the existing differences are due to this issue. In order to ensure the undeniable effect of the properties of these materials, modeling has been done for several samples, and their nonlinear analysis has been compared with the existing equation. In order to investigate the effect of different moduli of elasticity in these two analyses, materials are modeled in the software with a constant modulus of elasticity. Instead of three different moduli of elasticity, the longitudinal modulus has been used, which has a greater effect on the buckling. Additionally, this modulus of elasticity is used in Equation (1). Analyses pertaining to several samples are depicted in Figure 14 and Figure 15. The figures, labeled with the suffix “23,000” after the sample name, correspond to the results of the samples with a longitudinal modulus of 23,000 MPa.

Table 9. Results of isotropic modeling.

Sample Name	Buckling Load from Linear Analysis (N)	Buckling Load from Nonlinear Analysis (N)	Buckling Load from Equation (1) (N)	Difference between Linear Response and Equation (1) (%)	Difference between Nonlinear Response and Equation (1) (%)
IS1	3330.5	2633.46	2634.3	20.9	0.02
IS2	4896.7	4000	3488.3	28.8	13
IS3	4923.3	6785	6884.9	−39.8	−1
IS4	4929.3	8955	8304	−68	7.8
IS5	4934	6796	6019	−22	12.9
IS6	2279.3	1780	1748.8	23	1
IS7	3204	2644	2256.6	29.6	17
IS8	3256.4	4834	4397.8	−35	9.9
IS9	3267.1	5420	5026.9	−53.9	7.8
IS10	3278.7	3965	3722.9	13	6.5

presents the results of isotropic modeling for several samples. The ‘IS’ followed by a number designates the isotropic modeling for Sample ‘S’ with the corresponding number. According to this table, the difference between the nonlinear response and the response resulting from Equation (1) for 90% of the modeled samples is less than 13%. The difference between the response of nonlinear analysis and the equation is 0.02% in the best case and 17% in the worst case. Therefore, the use of an elasticity modulus alone in the equation includes an important part of the difference between the results of the analytical and numerical modeling. To correct this equation, a coefficient of transverse and height modulus of elasticity can also be included. Also, according to the graphs drawn in this section, it can be concluded that other parameters such as the ratio of length to restrained length, ratio of length to cross-sectional height, as well as thickness can also be factors affecting the moment and buckling load of the studied beams. Entering these parameters with appropriate coefficients can bring the response obtained from the equation and finite element analysis closer and can make the use of the equation wider for non-isotropic materials.

4.3. Key Findings and Implications

This study marks a significant departure from conventional research in structural engineering, which primarily focuses on beams without restraints. These findings demonstrate that the inclusion of lateral restraints in PFRP beams significantly enhances structural safety and stability. This discovery is pivotal, suggesting a necessary shift from traditional design approaches that may not fully capitalize on the stability benefits offered by lateral restraints. Moreover, this research indicates that PFRP beams, with their improved buckling performance due to lateral restraints, are prime candidates for wider adoption in environments prone to corrosion. A crucial aspect of this study is the recognition that conventional equations formulated for isotropic materials like steel are not entirely applicable to PFRP beams with restraints. Therefore, there is an evident need for the development of specific design standards that accurately reflect the unique characteristics of PFRP materials. These insights not only contribute to the academic understanding of PFRP beams but also have substantial implications for practical applications in structural engineering, heralding a new era of bridge design and construction.

5. Conclusions

In this study, the torsional buckling behavior of pultruded channel beams in bridges under pure bending is studied numerically. This research is among the first to analyze the impact of lateral restraints on PFRP beams, contributing significantly to the understanding and application of buckling behavior control in structural engineering. Non-isotropic behavior is chosen for the finite element modeling for the beams. The available analytical results have been compared with the results of finite element analysis for channel beams under pure bending. The differences between the values obtained from the software and the analytical results for PFRP channel beams indicate that the analytical equations used for these beams need to be modified to be able to more accurately estimate the buckling load of FRP beams under the conditions investigated in this study. It is worth mentioning that the scope of the results of this study is limited to the situations considered for analytical samples in the research, but it is expected that these results will have a more comprehensive application.

• The instability mode of all samples was caused by lateral-torsional buckling, and they started to twist and move laterally from the beginning of loading, and when the load reached the critical limit, they experienced lateral-torsional buckling. In fact, the buckling of all samples is gradual and accompanied by lateral deformation and twisting in the longitudinal direction. The reason for this issue was the existence of initial imperfections in PFRP beams due to nonlinear analysis, which affects the type of buckling.
• The application of restraint has resulted in the buckling load being magnified, with an augmentation factor ranging from at least 1.2 times to a maximum of 3.3 times that of the corresponding unrestrained sample. The smallest increase is observed when the restraint is positioned closest to the support, while the largest occurs when the restraint is situated in the middle of the span.
• Through a comparative examination of nonlinear analyses on samples sharing the same critical length yet possessing varying span lengths, it becomes evident that, as the beam span length approaches the critical length, the buckling load increases. Additionally, the enhancement in thickness corresponds to a diminishing effect of the span length differential in relation to the restrained span length on the buckling load.
• As the ratio of length to height in the section increases, the buckling load experiences a nonlinear reduction. Likewise, as thickness is enhanced, the influence of the increased ratio of length to cross-sectional depth on the buckling load diminishes.
• As the ratio of the length to the length of the restrained span increases, the buckling load exhibits a nonlinear increase, while concurrently, the effect of increased thickness on the critical buckling load decreases.
• In certain specimens, the outcome of nonlinear analysis exceeded the anticipated value, while in other cases, the result of buckling analysis surpassed expectations. This phenomenon can be attributed to the development of a form of rigidity and relative clamping at the end of the beam, induced by the presence of lateral restraints near the roller support. Consequently, this effect leads to an augmentation in the critical moment magnitude.
• The buckling load resulting from Equation (1), which is stated in the EN 1993 Eurocode 3 regulation [34], and the nonlinear response approach each other as the length-to-height ratio of the beam section is increased. Given the consistency of the cross-sectional height across all samples, it can be concluded that, as the length of the beam increases, the orthotropic behavior of the beam approaches the isotropic behavior of the beam and has less of an impact on the buckling load.
• The isotropic equation affords a closer approximation for beams with a greater restrained length. It can be asserted that the isotropic equation for orthotropic elastic beams yields a more accurate response, and the orthotropic behavior of beams exerts a more significant impact on the buckling load at lesser restraining distances.

6. Future Research

The significant gap between the results obtained from the equation and the nonlinear analysis underscores the limited accuracy of the employed equation in characterizing these materials. This highlights the pressing need for a refined equation that can offer a more precise estimation of the buckling load for the FRP beams under investigation in this research.

• Experimental validation: Conduct experimental studies to validate the numerical findings, particularly the impact of lateral restraints on the buckling behavior of PFRP beams.
• Extension to other beam types: Explore the effect of lateral restraints on other types of beams, such as I-section or tubular profiles, to broaden the understanding of buckling behavior in different structural elements.
• Long-term performance studies: Investigate the long-term performance of PFRP beams with lateral restraints under various environmental conditions, including temperature variations and exposure to moisture and chemicals.
• Development of improved analytical models: Develop and validate new analytical models that more accurately reflect the buckling behavior of PFRP beams, considering factors such as lateral restraints and material anisotropy.
• Integration with sustainable design practices: Incorporate findings into sustainable design practices, particularly for structures in corrosive environments or where material longevity is a concern.