Effect of bending test setup on the structural stability of I-shape pultruded FRP profiles using variable step size FE analysis

Abstract

Despite the widespread use of pultruded Fibre-reinforced Polymers (FRP) profiles, there is still a lack of standardised guidelines for the design of bending test setups on the structural-size level to match the structural stability demands of real applications in construction. This study investigates and compares the structural instability failure modes of lateral-torsional buckling, local buckling, and web crippling of I-shape pultruded FRP profiles under three-point and four-point bending. A novel modelling approach capable of capturing the longitudinal and sectional instabilities was established and verified against published experimental data and against the theory. The validated modelling approach was used to perform parametric studies on the bending test setup parameters, namely the span-to-depth ratio L/d and shear span-to-depth ratio a/d. The effect of these parameters on the failure modes and their thresholds was investigated and their interactions were characterised. Guidelines and recommendations were proposed for consideration when designing test setups and span lengths of I-shape pultruded FRP beams to widen the use of these profiles with more confidence in broader structural applications.

Keywords:

• I-shape pultruded GFRP beam
• lateral-torsional buckling
• local buckling
• web crippling
• bending test setup

1. Introduction

Pultruded Fibre-reinforced Polymers (FRP) profiles are growingly utilised in construction and infrastructure sectors due to their excellent corrosion resistance, cost-effective manufacturing process, and high strength-to-weight ratio (Aldoumani et al., 2016; Onyelowe et al., 2022). I-shape profiles represent more than 40.6% of the used pultruded FRP sections in civil structural applications, such as beams and panels in buildings and girders and decks in marine structures and bridges (Alhawamdeh, Alajarmeh, Aravinthan, Shelley, Schubel, Mohammed, et al., 2021). These profiles are prone to structural instability failure modes under compressive stresses, such as Euler buckling, lateral-torsional buckling, local buckling, and web crippling due to their application-driven slenderness and low anisotropic stiffness (F. Ascione et al., 2015, 2016). The bending behaviour, structural stability, and failure mode of pultruded FRP beams depend on the load arrangement, boundary conditions, and profile dimensions, all of which are significant parameters when setting up the test configuration (Alhawamdeh, Alajarmeh, Aravinthan, Shelley, Schubel, Mohammad, et al., 2021). The two most common bending test setups documented in the literature are the three-point and four-point bending shown in Figure 1. The span-to-depth ratio (L/d) in three-point bending and the shear span-to-depth ratio (a/d) along with the L/d ratio in four-point bending are key parameters responsible for configuring these test setups and controlling the results (Muttashar et al., 2016). In four-point bending, the shear span (a) is defined as the clear distance between the supporting point and the adjacent loading point at which the shear force remains constant (Xin et al., 2017).

Figure 1. Typical bending test setups: (a) three-point bending and (b) four-point bending.

However, there is still a lack of knowledge and standardised guidelines for the design of bending test setups on the structural-size level of these profiles. It is important to accurately characterise the structural behaviour and failure modes of these beams using the appropriate test setups that match the loading and supporting conditions and dimensions of the targeted applications in construction. Investigating the effect of the bending test setup parameters (L/d and a/d ratios) is essential to develop standardised provisions and guidelines for test setup configurations for these profiles on the structural-size level. This will widen the use of these profiles in broader civil structural applications with more confidence in their characteristic properties resulting from these tests.

The available standards and guidelines on bending test setup details present recommendations and configurations to characterise the flexural behaviour of FRP composites for specific applications and limited failure modes (ASTM D7249/D7249M–20, 2020; ASTM D8019–15, 2016). The ASTM D8019–15 (ASTM D8019–15, 2016) describes the test setup used to estimate the flexural strength of FRP cross-arms using three-point bending. The considered failure mode is local buckling of the compression flange with a span-to-depth ratio (L/d) > 16.0. The ASTM D7249/D7249M–20 (ASTM D7249/D7249M–20, 2020) reports the test setup used to assess the flexural strength of sandwich FRP beams and panels using four-point bending with a span-to-depth ratio (L/d) > 20.0. Following these standards, more studies and guidelines are needed to cover all the possible structural instability failure modes of structural-size profiles and report the effect of the associated test setup parameters on these failure modes comprehensively.

Many studies investigated the structural behaviour of pultruded FRP profiles subjected to three-point (Borowicz & Bank, 2011; Laudiero et al., 2013; Madenci et al., 2020) or four-point (Alhayek et al., 2022; Guades et al., 2014; Singh & Chawla, 2018) bending. However, only few studies examined the effect of the test setup parameters (L/d and a/d ratios) on the structural instability failure modes. The effect of the shear span-to-depth ratio was studied experimentally for hollow box pultruded FRP beams subjected to four-point bending (Muttashar et al., 2016). The investigated range was 1.2 ≤ a/d ≤8.0 with local buckling of the top flange and crushing of the thick web (web height/web thickness = 14.2) as the two dominant failure modes. The flexural behaviour of the experimentally tested beams was significantly affected by the shear span-to-depth ratio as increasing the a/d ratio (a/d > 7.5) increased the bending stress contribution by 85% while decreasing the a/d ratio (a/d < 5.0) increased the shear stress contribution by 50%. The flexural and torsional moduli of I-shape pultruded FRP profiles subjected to three-point and four-point bending were investigated experimentally for a range of span-to-depth ratios (Roberts & Al-Ubaidi, 2002). Three different cross-sectional sizes were studied with 7.5 ≤ L/d ≤30.0 for three-point bending and 19.0 ≤ L/d ≤30.0 and 3.2 ≤ a/d ≤7.5 for four-point bending. A close correlation was found between the flexural and torsional rigidities obtained from these test setups but the associated failure modes were not studied. The flexural behaviour of pultruded FRP sheet piles subjected to four-point bending was investigated for two different span-to-depth ratios of 8.0 and 17.0 (Wang et al., 2016). The influence of the L/d ratio was insignificant on the load capacity because of the localised debonding failure between the joined pairs. However, the deformation increased significantly for long-span specimens. The flexural capacity of hollow square and rectangular GFRP profiles was studied under three-point bending for various span-to-depth ratios (Estep et al., 2016). An analytical model was proposed to predict the failure load based on the L/d ratio. The accuracy of the model varied from high when 10.0 < L/d ≤42.0 with a bending failure in the tension flange to low when 4.0 ≤ L/d < 6.0 with a local buckling in the compression flange. Most of the available literature discussed the effect of L/d and a/d ratios on the test setups without covering the entire design space of these parameters. Moreover, the effect of these parameters on all the structural instability failure modes of I-shape FRP beams subjected to three-point and four-point bending was not reported holistically.

Conducting a parametric study to investigate the entire design space of the test setup parameters and their effect on the failure modes is not feasible using experimental or analytical approaches due to the limited experimental resources and unavailable closed-form equations, respectively. Consequently, a literature review was undertaken to explore the suitable numerical modelling approach for investigating the bending test setup parameters and their effect on the structural instability of I-shape pultruded FRP profiles. Numerical approaches such as the Finite Strip Method (FSM), Generalised Beam Theory (GBT), and Finite Element Method (FEM) are generally preferred (Kreja, 2011; Trang et al., 2022). FEM represents the suitable approach to simulate the bending test setups of I-shape pultruded FRP beams due to its superior flexibility and accuracy in modelling complex failure modes and versatile boundary and loading conditions (Alhawamdeh, Alajarmeh, Aravinthan, Shelley, Schubel, Mohammed, et al., 2021; Koloor SS et al., 2018; Lamberti et al., 2022).

However, commercial FEM software programs do not provide one analysis type convenient for modelling all types of structural instabilities (Xu et al., 2013). Alternatively, specific analysis methods are usually preferred depending on the targeted failure mode whether it is localised across the section or extended along the profile length (Sayyad & Ghugal, 2017). Moreover, the mechanism of handling the geometric imperfections during the FE analysis significantly affects the resulting buckling behaviour of the pultruded profile (F. Ascione, 2014). For example, the modified Riks (Arc-length) method is preferred for global instabilities such as global buckling and lateral-torsional buckling due to its load-incremental approach. Nevertheless, this method faces difficulties when simulating local instabilities like local buckling and web crippling because such sectional instabilities require the bifurcation point (i.e. buckling point) to be predefined by the user (Alhawamdeh, Alajarmeh, Aravinthan, Shelley, Schubel, Kemp, et al., 2021). The explicit dynamic analysis and Newton method (dynamic implicit analysis) are more suitable to track local instabilities such as local buckling and web crippling. However, the explicit dynamic analysis comprises some numerical parameters (e.g. artificial damping and mass scaling coefficients) that should be estimated correctly by experience to ensure accurate results (Hale & Ng, 2021). In addition, the dynamic implicit analysis requires the geometric imperfection history to locate the bifurcation point (L. Ascione et al., 2015). These limitations on different analysis methods raise the need for a comprehensive numerical approach capable of capturing both the global and local structural instabilities along the profile length and across the profile section. This numerical modelling approach is developed in this research based on the FEM with a variable step size algorithm.

This study aims to investigate the structural instability failure modes of I-shape pultruded FRP profiles subjected to three-point and four-point bending. A novel modelling approach was established to capture lateral-torsional buckling, local buckling, and web crippling using artificially imposed variable geometric imperfection functions guided by scaled shape functions capturing the longitudinal and sectional deformations and instabilities. The validated modelling approach was used to perform parametric studies on the bending test setup parameters, namely the span-to-depth ratio L/d and shear span-to-depth ratio a/d. The effect of these parameters on the structural instability failure modes was investigated and quantified. The thresholds between the instability failure modes were characterised graphically and statistically using failure maps and regression analysis. In addition, the interactions between the anisotropy ratio ($E_L/E_T$) and test setup parameters were studied. These parametric studies aimed to introduce the appropriate test setup configurations to match the structural stability demands of real applications in construction. These parametric studies provide guidelines and recommendations for accurate testing and characterisation of the bending behaviour and structural instability failure of I-shape pultruded FRP beams. This will widen the use of these profiles with more confidence in broader civil structural applications.

2. Numerical modelling approach

2.1. Material model and mesh

The material model including the elastic behaviour and the progressive failure behaviour was assembled using Abaqus/CAE 2019. The elastic behaviour of the FRP profile was simulated using elastic lamina material definition since it is compatible with laminated shells defined by two-dimensional plane stress formulation (Simulia, 2019). The input stiffness properties of the lamina consisted of the modulus of elasticity in the longitudinal direction (E₁), the modulus of elasticity in the transverse direction (E₂), the in-plane shear modulus of elasticity (G₁₂), the out-of-plane shear modulus of elasticity (G₂₃), and the in-plane Poisson's ratio ($v_{12}$). Since unidirectional laminae are considered transversely isotropic materials (Barbero, 2017), the value of the out-of-plane shear modulus (G₁₃) and (G₁₂) were assumed to be equal.

The progressive failure was simulated by a nonlinear analysis of geometry implementing the large displacement formulation and the Newton incremental method in the standard solver (Simulia, 2019) to capture all nonlinearities from the material failure, boundary conditions, and geometry. The Hashin damage model (Simulia, 2019) was used to define the progressive failure of fibres and matrix (resin). This model considers four failure modes: fibre failure in tension (rupture) and compression (buckling and kinking) and matrix failure in tension (cracking under transverse tension and shearing) and compression (crushing under transverse compression and shearing), as shown in Table 1. Each failure mode consists of a damage initiation criterion and a damage evolution (progressive failure) criterion.

Table 1. Failure modes and associated damage formulae utilised in the Hashin damage model Simulia, (2019)

Failure Mode	Damage Initiation^a	Damage Evolution^b
Fibre tension ($\mathit{\sigma }{\text{^}}_{11}$ ≥ 0)	${F_f}^{\mathit{t}}={\left(\frac{\mathit{\sigma }{\text{^}}_{11}}{X^T}\right)}^2+\mathit{\alpha }{\left(\frac{\mathit{\tau }{\text{^}}_{12}}{S^L}\right)}^2$	${{\mathit{\delta }}^{\mathit{ft}}}_{\mathit{eq}}=L^c\sqrt{{⟨{\mathit{\epsilon }}_{11}⟩}^2+\mathit{\alpha }{\left({\mathit{\epsilon }}_{12}\right)}^2},$ ${{\mathit{\sigma }}^{\mathit{ft}}}_{\mathit{eq}}=\frac{{\mathit{\sigma }}_{11}{\epsilon }_{11}+\mathit{\alpha }{\tau }_{12}{\epsilon }_{12}}{{{\mathit{\delta }}^{\mathit{ft}}}_{\mathit{eq}}/L^c}$
Fibre compression ($\mathit{\sigma }{\text{^}}_{11}$ < 0)	${F_f}^{\mathit{c}}={\left(\frac{\mathit{\sigma }{\text{^}}_{11}}{X^c}\right)}^2$	${{\mathit{\delta }}^{\mathit{fc}}}_{\mathit{eq}}=L^c⟨-{\epsilon }_{11}⟩,$ ${{\mathit{\sigma }}^{\mathit{fc}}}_{\mathit{eq}}=\frac{⟨-{\mathit{\sigma }}_{11}⟩⟨-{\epsilon }_{11}⟩}{{{\mathit{\delta }}^{\mathit{fc}}}_{\mathit{eq}}/L^c}$
Matrix tension ($\mathit{\sigma }{\text{^}}_{22}$ ≥ 0)	${F_m}^{\mathit{t}}={\left(\frac{\mathit{\sigma }{\text{^}}_{22}}{Y^T}\right)}^2+{\left(\frac{\mathit{\tau }{\text{^}}_{12}}{S^L}\right)}^2$	${{\mathit{\delta }}^{\mathit{mt}}}_{\mathit{eq}}=L^c\sqrt{{⟨{\mathit{\epsilon }}_{22}⟩}^2+{\left({\mathit{\epsilon }}_{12}\right)}^2},$ ${{\mathit{\sigma }}^{\mathit{mt}}}_{\mathit{eq}}=\frac{⟨{\mathit{\sigma }}_{22}⟩⟨{\epsilon }_{22}⟩+{\tau }_{12}{\epsilon }_{12}}{{{\mathit{\delta }}^{\mathit{mt}}}_{\mathit{eq}}/L^c}$
Matrix compression ($\mathit{\sigma }{\text{^}}_{22}$ < 0)	${F_m}^{\mathit{c}}={\left(\frac{\mathit{\sigma }{\text{^}}_{22}}{2S^T}\right)}^2+\left({\left(\frac{Y^C}{2S^T}\right)}^2-1\right)\frac{\mathit{\sigma }{\text{^}}_{22}}{Y^C}+{\left(\frac{\mathit{\tau }{\text{^}}_{12}}{S^L}\right)}^2$	${{\mathit{\delta }}^{\mathit{mc}}}_{\mathit{eq}}=L^c\sqrt{{⟨-{\mathit{\epsilon }}_{22}⟩}^2+{\left({\mathit{\epsilon }}_{12}\right)}^2},$ ${{\mathit{\sigma }}^{\mathit{mc}}}_{\mathit{eq}}=\frac{⟨-{\mathit{\sigma }}_{22}⟩⟨-{\epsilon }_{22}⟩+{\tau }_{12}{\epsilon }_{12}}{{{\mathit{\delta }}^{\mathit{mc}}}_{\mathit{eq}}/L^c}$

^a${F_f}^{\mathit{t}}$, ${F_f}^{\mathit{c}}$, ${F_m}^{\mathit{t}}$, and ${F_m}^{\mathit{c}}$ are damage initiation factors for fibre under tension and compression and matrix under tension and compression, respectively. α is the shear stress counter affecting the fibre tensile initiation criterion. $\mathit{\sigma }{\text{^}}_{11}$, $\mathit{\sigma }{\text{^}}_{22}$, and $\mathit{\tau }{\text{^}}_{12}$ refer to the axial, transverse, and shear stress components, respectively, of the effective stress tensor ${\sigma }^{\text{^}}$ responsible for evaluating the initiation criteria. ^b ${{\mathit{\delta }}^{\mathit{ft}}}_{\mathit{eq}}$, ${{\mathit{\sigma }}^{\mathit{ft}}}_{\mathit{eq}}$, ${{\mathit{\delta }}^{\mathit{fc}}}_{\mathit{eq}}$, ${{\mathit{\sigma }}^{\mathit{fc}}}_{\mathit{eq}}$, ${{\mathit{\delta }}^{\mathit{mt}}}_{\mathit{eq}}$, ${{\mathit{\sigma }}^{\mathit{mt}}}_{\mathit{eq}}$, ${{\mathit{\delta }}^{\mathit{mc}}}_{\mathit{eq}}$, and ${{\mathit{\sigma }}^{\mathit{mc}}}_{\mathit{eq}}$ are the equivalent displacement and stress components of the element under fibre tension and compression and matrix tension and compression, respectively. $L^c$ is the characteristic length of the element.

The strength limits controlling the damage initiation criteria comprise the longitudinal compressive ($X^c$) and tensile ($X^T$) strengths, the transverse compressive ($Y^C$) and tensile ($Y^T$) strengths, and the shear strengths in the longitudinal ($S^L$) and transverse ($S^T$) directions. When the damage initiation criterion is met for a specific failure mode (the related factor reaches a value of 1.0), the damage evolution criterion is triggered for the related failure mode. The damage evolution algorithm calculates the area under the stress-displacement curve of the element, which is equivalent to the energy release. The fracture energy components of the longitudinal compressive ($G_{\mathit{LC}}$) and tensile ($G_{\mathit{LT}}$) fracture energies and the transverse compressive ($G_{\mathit{TC}}$) and tensile ($G_{\mathit{TT}}$) fracture energies determine the ultimate energy required to cause the collapse of the element.

The damage stabilisation scheme in Abaqus/CAE 2019 was applied along with the Hashin damage model to convert the tangent stiffness matrix into a positive equivalent for sufficient time steps to prevent analysis difficulties, divergence, and errors accompanying the modelling of material failure. The damage evolution algorithms are damped by relaxing the time in the vicious system through a coefficient of viscosity introduced by the solver (Simulia, 2019).

The FRP profile was simulated using the in-plane general-purpose quadrilateral (8-node) continuum shell (SC8R) element. This element is rendered as a three-dimensional geometry allowing for a more accurate through-thickness seeding, modification, and analysis compared to conventional shells (Simulia, 2019). Thus, it represents one of the best options to model laminated shells. The hourglass modes (i.e. inaccurate deformation of the mesh caused by artificial energy generated in the structure) were mitigated and the kinematic response was accurately captured by increasing the number of elements through the thickness. The optimal element size to achieve accurate and converged results with the available computational resources was determined through a mesh sensitivity study. An element size with an edge length of 4 mm was selected along with four elements through the thickness.

2.2. Continuously variable step size for structural instability

Generally, geometric imperfections are used with multiple roles in FEM applications. One of their applications is to simulate the effect of out-of-straightness and manufacturing imperfections on the structural stability of FRP profiles by introducing the imperfections, measured experimentally or provided by the manufacturer, as initial nodal displacements in the model (Nguyen et al., 2015). Geometric imperfections can also be used for pure numerical purposes to assist the FE solver in locating the buckling point (Nguyen et al., 2023). Such imperfections are usually extracted as buckling modes from eigenvalue buckling analysis and then used as initial nodal displacements in the first step of a following non-linear buckling analysis to simulate a specific buckling type (Nguyen et al., 2013). In these applications, the geometric imperfections are defined and scaled by the user before the analysis to target a specific buckling type. In addition, these imperfections are introduced only once in the first step of non-linear buckling analysis with no further assessment or update during the analysis. Alternatively, the new modelling approach was designed to track the structural instability failure modes of lateral-torsional buckling, local buckling, and web crippling using variable geometric imperfection functions guided by scaled shape functions capturing the longitudinal and sectional deformations and instabilities. The modelling approach overrides the step size of time calculated by the FE solver and replaces it with a Continuously Variable Step Size (CVSS) at the beginning of each step. The CVSS was implemented and automated as a code in the FEM script using Python 3.9.1 programming language due to the compatibility between this language and the Abaqus kernel (i.e. the script interpreter) (Alhawamdeh et al., 2022). Abaqus/standard was used as an FE solver to undertake the nonlinear buckling analysis as discussed in the previous section.

The CVSS is calculated based on an artificially imposed geometric imperfection function ѱ(x,y,z) that considers the effect of the profile dimensions and bending test setup parameters (L/d and a/d ratios). The ѱ(x,y,z) function consists of three shape functions tracking the deformation of the beam across the section and along the span. ɸ_b(x,y,z), ɸ_f(x,y,z), and ɸ_w(x,y,z) represent the shape functions defined to trace the deformation of the beam, flange, and web neutral axes, respectively, as shown in Figure 2. Each shape function is defined as a path of nodes in the FE model and is assigned to capture a specific failure mode. The ɸ_b(x,y,z), ɸ_f(x,y,z), and ɸ_w(x,y,z) functions are designed to capture the lateral-torsional buckling of the beam, local buckling of the top flange, and web crippling, respectively Appendix A.

Figure 2. Three shape functions capturing the structural instability of the I-shape pultruded FRP beam.

The geometric imperfection function of the FRP profile is represented mathematically as the summation of the shape functions scaled depending on the intensity of their relevant deformations:

(1) $\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{b}\left(\mathit{t}\right)\mathit{b}\left(\mathit{x},\mathit{y},\mathit{z}\right)+\mathit{f}\left(\mathit{t}\right)\mathit{f}\left(\mathit{x},\mathit{y},\mathit{z}\right)+\mathit{w}\left(\mathit{t}\right)\mathit{w}\left(\mathit{x},\mathit{y},\mathit{z}\right)$(1)

where $\left(\mathit{x},\mathit{y},\mathit{z}\right)$ is the geometric imperfection function. ɸ_b(x,y,z), ɸ_f(x,y,z), and ɸ_w(x,y,z) are the shape functions of beam, flange, and web, respectively. ɣ_b(t), ɣ_f(t), and ɣ_w(t) are the shape function weights of the beam, flange, and web, respectively, as functions of the analysis time. These weights are updated from the path of field nodes at the beginning of each step and are calculated as follows:

(2) ${\mathit{\gamma }}_{\mathrm{b}}\left(\mathrm{t}\right)=\frac{{\mathit{\delta }}_{\mathrm{b}}}{{\mathit{\delta }}_{\mathrm{b}}+{\mathit{\delta }}_{\mathrm{f}}+{\mathit{\delta }}_{\mathrm{w}}}$(2)

(3) ${\mathit{\gamma }}_{\mathrm{f}}\left(\mathrm{t}\right)=\frac{{\mathit{\delta }}_{\mathrm{f}}}{{\mathit{\delta }}_{\mathrm{b}}+{\mathit{\delta }}_{\mathrm{f}}+{\mathit{\delta }}_{\mathrm{w}}}$(3)

(4) ${\mathit{\gamma }}_{\mathrm{w}}\left(\mathrm{t}\right)=\frac{{\mathit{\delta }}_{\mathrm{w}}}{{\mathit{\delta }}_{\mathrm{b}}+{\mathit{\delta }}_{\mathrm{f}}+{\mathit{\delta }}_{\mathrm{w}}}$(4)

Where δ_b, δ_f, and δ_w represent the normalised displacement of the beam, flange, and web shape functions, respectively. The normalised displacement is obtained by dividing the shape function displacements by the maximum displacement attained by the shape function at time t. This displacement represents a quantitative measure of the deformation intensity of each shape function. Contrary to the default analysis methods available in commercial FEM software, the proposed modelling approach guides the solver to the bifurcation point (i.e. buckling point) using the shape functions and their weights. The default analysis methods are eye-blinded and need the user’s assistance with providing predefined constant geometric imperfections to lead the solver to the accurate failure mode. Conversely, the proposed approach provides an automated tool to compute the variation in the geometric imperfections along the analysis time and update the step size by tracking the deformation of the shape functions and calculating their weights. If lateral-torsional buckling is the dominant failure mode, the normalised displacement of the beam shape function (δ_b) becomes large compared to the normalised displacements of the flange (δ_f) and web (δ_w). Consequently, the weight of the beam shape function (ɣ_b(t)) increases compared to the weights of the flange (ɣ_f(t)) and web (ɣ_w(t)) shape functions. Thus, the geometric imperfection function ѱ(x,y,z) will be updated with a larger out-of-straightness located along the beam length compared to a lower out-of-straightness across the beam section. This will obtain a larger (faster) step size after each step as the solution converges towards the accurate failure mode. The same procedure applies if the local buckling of the flange or web crippling is the dominant failure mode. Figure 3 depicts the typical response of the shape function weights versus the analysis time when each structural instability failure mode (lateral-torsional buckling, local buckling, or web crippling) dominates the solution.

Figure 3. Typical response of shape function weights versus the analysis time: (a) when lateral-torsional buckling is the dominant failure mode, (b) when local buckling is the dominant failure mode, and (c) when web crippling is the dominant failure mode.

The step size is adjusted after each step based on the updated geometric imperfection function, which guides the solver closer to the buckling point. When the updated geometric imperfection function attains a closer shape to the accurate failure mode, the analysis converges smoothly towards the solution with a high weight of that specific failure mode. Thus, the solver increases the step size to speed up the analysis with high certainty. Conversely, when the analysis diverges slowly from the solution, the weight of the shape functions is decreased. Consequently, the solver decreases the step size to slow down the analysis as the certainty of the results is low. This will enable the solver to accurately locate the convergence trend (i.e. locate the buckling point) benefiting from the small time steps. When the solution convergence is reattained, the weight of the relevant shape function increases again causing the solver to increase the step size again (speed up the analysis). Figure 4 summarises the new modelling approach introduced in this study incorporating the CVSS algorithm and Abaqus finite element analysis. This modelling scheme was combined and automated using Python 3.9.1 programming language because of its compatibility with Abaqus kernel (i.e. script interpreter) (Alhawamdeh et al., 2022).

Figure 4. Modelling scheme of variable step size FEM analysis.

This new modelling approach has two main advantages as it (i) efficiently locates the buckling (bifurcation) point by introducing updated geometric imperfections to guide the solver to the correct buckling type depending on the longitudinal and transverse deformations of the beam and (ii) reduces the computational cost (analysis time) by continuously adjusting the step size according to the solution convergence to avoid unnecessary small time steps. The proposed approach can simulate the structural instability failure modes of lateral-torsional buckling, local buckling, and web crippling under different bending test setup configurations. In the following section, three case studies of experimental tests were selected from the literature to verify the validity of the proposed approach.

3. Validation case studies

3.1. Background

Three experimental programs, each one addressing one structural instability failure mode, from published literature were considered as case studies to verify the proposed modelling approach. The first experimental study investigated the lateral-torsional buckling of I-shape pultruded FRP beams and was performed by (Nguyen et al., 2014). The second and third experimental studies examined the local buckling and web crippling of these beams and were conducted by (Yang et al., 2020) and (Fernandes et al., 2015), respectively. Table 2 summarises the dimensions of the specimens and the test setup details of the three validation case studies.

Table 2. Dimensions and test setup details of the validation case studies

Specimen ID^a	Cross-sectional dimensions (mm)^b	Beam length (mm)	Type of test setup	Reference
I-LTB	120.0×60.0×6.0	4064	Three-point bending	Nguyen et al. (2014)
I-LB	152.4×101.6×6.35	2900	Four-point bending	Yang et al. (2020)
I-WC	120.0×60.0×6.0	240	Three-point bending	Fernandes et al. (2015)

^aI-LTB refers to an I-shape pultruded FRP beam under lateral-torsional buckling. I-LB denotes an I-shape pultruded FRP beam under local buckling. I-WC points to an I-shape pultruded FRP beam under web crippling. ^b the cross-sectional dimensions of the I-shape beam were labelled as follows: web height × flange width × wall thickness.

The same boundary and loading conditions of the experimental tests were simulated in the FE models. Cylindrical solid parts with three-dimensional geometry, discrete rigid definition, and a radius of 20 mm were used to model the supporting and loading fixtures. For the boundary conditions, all translational degrees of freedom were fixed at the profile ends. Unrestrained rotation was permitted for the rotational degrees of freedom. Regarding the loading condition, a quasi-static loading rate of 10.0 mm/min was applied at the loading points using a vertical downward translational degree of freedom. The boundary and loading conditions were assigned on the cylindrical fixtures and were kinematically connected to the profile using a tie contact definition.

The mechanical properties of the pultruded FRP profiles tested in the case studies were extracted from the relevant references, as shown in Table 3. E_L and E_T refer to the longitudinal and transverse elastic moduli, respectively. G_LT and $v_{\mathit{LT}}$ denote the in-plane shear modulus and in-plane Poisson’s ratio, respectively. These values were used in the FEM analysis to simulate the structural behaviour of the I-LTB, I-LB, and I-WC I-shape beam specimens subjected to lateral-torsional buckling, local buckling, and web crippling failure modes, respectively. The intended experimental failure modes were guaranteed to be obtained as the same experimental test setup parameters (i.e. L/d and a/d) of the case studies were used in the FE models. All the profiles in these case studies were characterised by a high slenderness ratio of the web (web height/web thickness ≥20.0), which alleviated the shear deformation effects that usually occur in thick cross-sectional elements subjected to transverse loading, such as the web (Correia et al., 2011).

Table 3. Characteristic mechanical properties of the I-shape pultruded FRP profiles used in the validation case studies

Reference	Nguyen et al. (2014)	Yang et al. (2020)	Fernandes et al. (2015)
EL (MPa)	23,000	24,500	26,700
ET (MPa)	11,200	9,400	9,700
GLT (MPa)	3,000	2,900	3,900
$v_{\mathit{LT}}$	0.28	0.32	0.30

3.2. Results comparison

The validity of the new modelling approach was assessed by comparing the numerical results against the experimental results in terms of the failure mode and buckling strength. Moreover, the proposed CVSS analysis was compared with the default Abaqus/standard analysis in terms of solution convergence and analysis speed. Figure 5 depicts the experimental failure modes (cited from the aforementioned references) versus the CVSS FEM failure modes for all the considered structural instabilities of lateral-torsional buckling, local buckling of the flange, and web crippling. The failure modes resulted from the new modelling approach matched the experimental failure modes reported in the literature.

Figure 5. Experimental versus CVSS-FEM failure modes of the validation case studies: (a) lateral-torsional buckling of I-LTB (experimental photo cited from Nguyen et al. (2014)), (b) web crippling of I-WC (experimental photo cited from Fernandes et al. (2015)), and (c) local buckling of I-LB (experimental photo cited from Yang et al. (2020)).

The numerical buckling strength obtained from the CVSS FEM analysis was also compared against the experimental buckling strength extracted from the relevant literature. In addition, the theoretical buckling strength of each failure mode was calculated from the available closed-form equations (Bank, 2006) and was compared against the FEM results. Table 4 presents a comparison between these results of the I-shape FRP beam specimens for the same test setup parameters provided in the case studies. The percentile error between the FEM and the experimental and theoretical buckling strengths was 5.95% and 4.04%, respectively, for the specimen dominated by lateral-torsional buckling (I-LTB). Under local buckling failure, the error between the FEM and the experimental and theoretical buckling strengths of specimen I-LB was 6.47% and 4.33%, respectively. The percentile error between the FEM and the experimental and theoretical buckling strengths was 5.39% and 4.49%, respectively, for the I-WC specimen which failed by web crippling. These results emphasised the satisfactory agreement between the proposed modelling approach and the published data and theory.

Table 4. Numerical versus experimental and theoretical buckling strength results of the validation case studies of I-shape FRP beams subjected to bending

Specimen ID	Test setup parameters	Theoretical buckling strength (MPa)^a (Bank, 2006)	Experimental buckling strength (MPa)	FEM buckling strength (MPa)
I-LTB	L/d = 33.33 a/d = N/A	${\sigma }_{cr}^{lat}=\frac{C_b}{S_x}\sqrt{\frac{{\pi }^2{E}_L{I}_y{G}_{\mathit{LT}}\mathit{J}}{{\left(k_f{L}_b\right)}^2}+\frac{{\pi }^4{E}_L^2{I}_y{C}_w}{{\left(k_f{L}_b\right)}^2{\left(k_w{L}_b\right)}^2}}=47.82$	46.87 Nguyen et al. (2014)	49.84
I-LB	L/d = 19.03 a/d = 6.56	${\sigma }_{cr}^{\mathrm{local}}=\frac{1}{{\left(\frac{b_f}{2}\right)}^2{t}_f}\left(7\sqrt{\frac{D_L{D}_T}{1+4.12{\zeta }_{\mathit{flange}}}}+12D_S\right)=173.98$	170.08 Yang et al. (2020)	181.85
I-WC	L/d = 2.00 a/d = N/A	${\sigma }_{cr}^{\mathrm{crippling}}=\frac{2{\pi }^2}{t_w{b}_{eff}^2}\left(\sqrt{D_L{D}_T}+D_{\mathit{LT}}+2D_S\right)=223.35$	221.24 Fernandes et al. (2015)	233.85

^a${\sigma }_{cr}^{lat}$ represents the critical lateral-torsional buckling stress. $C_b$ and $C_w$ are coefficients for moment variation and warping along the beam, respectively. $\mathit{J}$, $I_y$, and $S_x$ are the torsional constant, moment of inertia about the y-axis, and section modulus, respectively. $k_f$ and $k_w$ are the effective length factor for flexural and torsional buckling, respectively. $L_b$ is the unbraced length of the beam. $E_L$ and $G_{\mathit{LT}}$ are the longitudinal elastic modulus and in-plane shear modulus, respectively. ${\sigma }_{cr}^{\mathrm{local}}$ is the critical local buckling stress. $b_f$ and $t_f$ are the flange width and thickness, respectively. $D_L$, $D_T$, and $D_S$ are the longitudinal, transverse, and shear flexural stiffness components of the laminated profile. ${\zeta }_{\mathit{flange}}$ is the coefficient of restraint between the flange and web. ${\sigma }_{cr}^{\mathrm{crippling}}$ represents the critical web crippling stress. $t_w$ and $b_{\mathit{eff}}$ are the web thickness and depth, respectively.

)

The CVSS FEM analysis was compared against the default Abaqus/standard analysis to examine the solution convergence and analysis speed, as shown in Figure 6. Both methods used the large displacement formulation and the Newton incremental method along with the elastic lamina material definition and Hashin damage model. The only difference between the two methods is the step size calculations at the beginning of each step, as discussed in section 2.2. Unlike the default analysis method, the CVSS analysis continuously examines and updates the step size based on a variable geometric imperfection function guided by scaled shape functions capturing the longitudinal and sectional deformations and instabilities.

Figure 6. Solution convergence of the CVSS analysis against default Abaqus/Standard analysis for the validation case studies.

When the number of increments (steps) in the analysis is increased, each step size is decreased resulting in a slow analysis with high computational cost. Conversely, when the number of increments is decreased, each step size is increased resulting in a fast analysis with low computational cost. The CVSS analysis was 1.24, 1.27, and 1.33 faster than the default analysis for the I-LTB, I-LB, and I-WC specimens dominated by lateral-torsional buckling, local buckling, and web crippling, respectively. The web crippling analysis of the I-WC specimen was the fastest followed by the I-LB specimen and then the I-LTB specimen. This was referred to as the effect of the model size (i.e. the number of elements in the model) on the computational efficiency since the size of the I-LTB FEM model was the largest followed by the I-LB and then the I-WC model. In terms of the buckling strength, it is noted that the new modelling approach obtains closer values to the experimental results with lower percentile error compared to the default FE analysis, as presented in Figure 6. This was referred to as the effect of the CVSS algorithm in simulating the geometric imperfections required for faster and more accurate buckling analysis.

This new modelling approach proved its capability and accuracy as a robust numerical tool that can capture the different structural instability failure modes present in bending tests. It can be used to perform parametric studies on the bending test setup parameters (i.e. span-to-depth ratio L/d and shear span-to-depth ratio a/d). In the following section, comprehensive parametric studies were undertaken on these parameters in order to conclude testing guidelines and recommendations for accurate characterisation of the bending behaviour and structural instability failure of pultruded FRP profiles. Consequently, the reliability of the characteristic properties of these profiles obtained based on the appropriate test setup configurations matching the demands of the intended applications will increase. This will widen the use of these profiles with more confidence in broader civil structural applications. It is worth noting that FRP profiles with variable mechanical properties across the section are out of scope in the following parametric studies. The proposed modelling approach can be applied to such profiles but with minor changes to the composite layup definition in Abaqus (i.e. defining each structural element in the cross-section as an independent part in the model) in order to assign the provided mechanical properties for each part of the cross-section accordingly.

4. Parametric study on the test setup parameters

In this section, the validated modelling approach was used to undertake parametric studies on the bending test setup parameters, namely the span-to-depth ratio L/d and shear span-to-depth ratio a/d. The effect of these parameters on the structural instability failure modes of lateral-torsional buckling, local buckling, and web crippling was investigated and quantified. The threshold (i.e. the inflection point at which the failure mode changes to another one) between the failure modes was characterised graphically and statistically using failure maps and regression analysis, respectively, to establish the relationship between the test setup parameters and the shift of failure mode. In addition, the interaction between the anisotropy of I-shape pultruded FRP profiles and test setup parameters was studied. These parametric studies were aimed to provide testing guidelines and recommendations for accurate characterisation of the bending behaviour and structural instability failure of pultruded FRP profiles.

4.1. Effect of span-to-depth ratio (L/d)

The results of the parametric study on the effect of the L/d ratio on the buckling capacity and failure mode of I-shape pultruded FRP beams were highlighted by generating failure maps. These failure maps are shown in Figures 7 and 8 for three-point and four-point bending tests, respectively. The buckling strength of each failure mode was normalised to provide universal representative results applicable to other I-shape profiles with different material strengths and layup properties. The normalised strength ($\mathit{\lambda }$) was introduced by dividing the buckling strength ${\sigma }_{\mathit{BL}}$ (obtained from a structural-level test) by the ultimate material strength ${\sigma }_u$ (obtained from a coupon test):

Figure 7. Failure map showing the effect of span-to-depth ratio on I-shape pultruded FRP profiles subjected to three-point bending with different anisotropy ratios.

Figure 8. Failure map showing the effect of span-to-depth ratio on I-shape pultruded FRP profiles subjected to four-point bending with different anisotropy ratios ($\mathit{a}/\mathit{d}=\frac{1}{4}\mathit{L}/\mathit{d}$).

(5) $\mathit{\lambda }={\sigma }_{\mathit{BL}}/{\sigma }_u$(5)

The ultimate material strength was considered as the average of the characteristic compressive strength values reported in the three validation case studies (342.0, 296.0, and 318.0 MPa from (Nguyen et al., 2014; Yang et al., 2020), and (Fernandes et al., 2015), respectively), which equals 318.7 MPa. The characteristic compressive strength was selected as the normalisation standard value since it controls the ultimate capacity of FRP laminated beams, which generally collapse by a compressive failure of fibres of the top flange under bending (Alhawamdeh, Alajarmeh, Aravinthan, Shelley, Schubel, Mohammad, et al., 2021). The generated failure maps were also checked against the published experimental data of the three validation case studies. These failure maps can be used to estimate the buckling strength and failure mode of I-shape FRP beams with different configurations of the bending test setup parameters without the need for large experimental programs but only a limited number of experiments to verify and correlate the results against these charts. The data presented in Figures 7 and 8 were used to statistically estimate the relationship between the span-to-depth ratio and anisotropy ratio and the normalised strength by performing a polynomial regression analysis using the statistical analysis software Minitab 19 (Minitab-19, 2017). The resulted regression formulae for three-point and four-point bending tests are presented in equations 6(6) $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.7964-0.0002\mathit{L}/\mathit{d}\pm 0.03E_L/E_T,\mathit{L}/\mathit{d}\le 11.0\\ 1.2398+0.0007{(\mathit{L}/d)}^2-0.0534\mathit{L}/\mathit{d}\pm 0.033E_L/E_T,11.0\le \mathit{L}/\mathit{d}\le 26.5\\ 2.1221+0.0014{(\mathit{L}/d)}^2-0.1056\mathit{L}/\mathit{d}\pm 0.021E_L/E_T,\mathit{L}/\mathit{d}\ge 26.5\end{array}\right.$(6) and 7(7) $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.8354-0.0009\mathit{L}/\mathit{d}\pm 0.032E_L/E_T,\mathit{L}/\mathit{d}\le 12.5\\ 2.1646+0.0034{(\mathit{L}/d)}^2-0.1544\mathit{L}/\mathit{d}\pm 0.028E_L/E_T,12.5\le \mathit{L}/\mathit{d}\le 24.0\\ 1.4902+0.001{(\mathit{L}/d)}^2-0.0713\mathit{L}/\mathit{d}\pm 0.017E_L/E_T,\mathit{L}/\mathit{d}\ge 24.0\end{array}\right.$(7) , respectively. Each formula consists of three polynomials each one describes a specific failure zone starting from lateral-torsional buckling, local buckling, and web crippling bounded by $\mathit{L}/\mathit{d}\ge 26.5$, $11.0\le \mathit{L}/\mathit{d}\le 26.5$, and $\mathit{L}/\mathit{d}\le 11.0$ for three-point bending and $\mathit{L}/\mathit{d}\ge 24.0$, $12.5\le \mathit{L}/\mathit{d}\le 24.0$, and $\mathit{L}/\mathit{d}\le 12.5$ for four-point bending (when $\mathit{a}/\mathit{d}=\frac{1}{4}\mathit{L}/\mathit{d}$), respectively. The coefficient of determination R² that describes the accuracy of the regression model was 0.93 and 0.91 for three-point and four-point bending, respectively.

(6) $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.7964-0.0002\mathit{L}/\mathit{d}\pm 0.03E_L/E_T,\mathit{L}/\mathit{d}\le 11.0\\ 1.2398+0.0007{(\mathit{L}/d)}^2-0.0534\mathit{L}/\mathit{d}\pm 0.033E_L/E_T,11.0\le \mathit{L}/\mathit{d}\le 26.5\\ 2.1221+0.0014{(\mathit{L}/d)}^2-0.1056\mathit{L}/\mathit{d}\pm 0.021E_L/E_T,\mathit{L}/\mathit{d}\ge 26.5\end{array}\right.$(6)

(7) $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.8354-0.0009\mathit{L}/\mathit{d}\pm 0.032E_L/E_T,\mathit{L}/\mathit{d}\le 12.5\\ 2.1646+0.0034{(\mathit{L}/d)}^2-0.1544\mathit{L}/\mathit{d}\pm 0.028E_L/E_T,12.5\le \mathit{L}/\mathit{d}\le 24.0\\ 1.4902+0.001{(\mathit{L}/d)}^2-0.0713\mathit{L}/\mathit{d}\pm 0.017E_L/E_T,\mathit{L}/\mathit{d}\ge 24.0\end{array}\right.$(7)

Unlike the lateral-torsional buckling and local buckling failure zones, the web crippling failure zone was characterised by a linear $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)$ function with a very small slope (approximately constant). It was inferred that the effect of the L/d ratio on the web crippling capacity is negligible considering the a/d ratio is constant. This was referred to the dependency of web crippling on the distance between the supporting and loading points (which is defined by the a/d ratio), while the lateral-torsional buckling and local buckling significantly depend on the L/d ratio as evidenced by the steep $\mathit{\lambda }\left(\mathit{L}/\mathit{d},E_L/E_T\right)$ curves in their failure zones.

Moreover, the threshold boundaries between the failure zones were shifted when comparing the effect of the L/d ratio on three-point versus four-point bending, as shown in Figures 7 and 8. The web crippling extended on a larger failure zone by 1.14 in four-point compared to three-point bending. This was referred to the smaller distance between the loading and supporting points in four-point bending compared to three-point bending for the same L/d ratio, which increased the transverse concentrated load on the web at the cost of decreasing the compressive stress on the top flange. Thus, the web crippling failure was extended over a larger zone and the local buckling failure zone was reduced under four-point bending. This behaviour led to a higher web crippling strength under four-point up to 3.9% compared to three-point bending as the load is distributed over two loading points in four-point bending, which alleviated the concentrated stress, compared to one loading point in three-point bending.

On the contrary, lateral-torsional buckling and local buckling extended on smaller failure zones by 0.91 and 0.88, respectively, in four-point compared to three-point bending. This behaviour in the lateral-torsional buckling zone was referred to the longer unbraced span (i.e. the distance between the supporting point and the nearest loading point) in three-point (that is L/2) compared to four-point bending (that is a) for the same L/d ratio. Consequently, the lateral-torsional buckling capacity was higher under four-point bending by up to 12.6% compared to three-point bending due to the two bracing points (loading points) in four-point bending against only one bracing point at mid-span in three-point bending. The smaller failure zone of local buckling in four-point bending was referred to the higher compressive stress concentrated on the top flange at one loading point in three-point bending compared to the lower compressive stress distributed on two loading points in four-point bending. The concentrated compressive stress on the top flange was higher by up to 7.2% (depending on the anisotropy ratio) under three-point bending, thus, the local buckling strength was higher under four-point bending by up to 7.2%.

Overall, the web crippling capacity was the highest followed by the local buckling capacity and then the lateral-torsional buckling capacity. This was referred to the inverse relationship between the nominal capacity of the beam and the span length. Another factor was the difference in the stress nature and distribution between these failure modes. In web crippling, the stress of concern across the section is bearing stress caused by the transverse compression loading. Conversely, the stress of concern in local buckling and lateral-torsional buckling is bending stress (longitudinal compressive stress on the top flange), with a torsional stress effect added in lateral-torsional buckling.

The effect of the laminate anisotropy ratio $E_L/E_T$ (i.e. the longitudinal-to-transverse elastic moduli ratio) was also examined in the failure maps, as shown in Figures 7 and 8. The lateral-torsional buckling strength increased up to 16.1% and 11.9% when the $E_L/E_T$ ratio increased from 1.0 up to 3.0 under three-point and four-point bending, respectively. Conversely, the local buckling and web crippling strengths increased up to 35.3% and 41.7% under three-point bending and up to 30.6% and 46.2% under four-point bending, respectively, when the $E_L/E_T$ ratio decreased from 3.0 up to 1.0. This was referred to as the effect of increasing $E_L$ on enhancing the longitudinal rigidity of the beam, which reduces the lateral displacement accompanying lateral-torsional buckling of the beam about its neutral axis. Thus, the lateral-torsional buckling strength is increased. On the other hand, increasing $E_T$ enhances the sectional rigidity of the flange and web, which reduces the sectional displacement accompanying local buckling of the flange or web crippling. Consequently, the local buckling or web crippling strength is increased. In addition, the effect of laminate anisotropy was decreased when the L/d ratio increased to reach the lateral-torsional buckling zone. This was referred to as the dominant influence of the beam span length on the lateral-torsional buckling strength compared to the laminate layup properties.

Moreover, the anisotropy ratio influenced the threshold boundaries between the failure zones. Increasing the $E_L/E_T$ ratio decreased the size of web crippling and local buckling failure zones up to 0.68 and 0.74 under three-point bending and up to 0.66 and 0.75 under four-point bending, respectively. This was referred to as the effect of the low value of transverse modulus (when increasing the $E_L/E_T$ ratio) on decreasing the sectional stability and the related buckling capacities. Thus, a beam with a low value of transverse modulus requires a shorter span to attain the same web crippling or local buckling capacity of a beam with a high value of transverse modulus and a longer span. On the contrary, decreasing the $E_L/E_T$ ratio enhanced the sectional stability and increased the web crippling and local buckling capacities. Consequently, the size of the related failure zones increased by up to 1.47 and 1.35 under three-point bending and up to 1.51 and 1.33 under four-point bending, respectively. The beam with higher values of transverse modulus and L/d ratio can attain the same buckling capacity of a beam with lower values of transverse modulus and L/d ratio.

4.2. Effect of shear span-to-depth ratio (a/d) on four-point bending

The effect of the shear span-to-depth ratio (a/d) on the normalised buckling strength ($\mathit{\lambda }$) and failure mode of I-shape beams subjected to four-point bending was examined for different anisotropy ratios, as shown in Figure 9. The statistical analysis software Minitab 19 (Minitab-19, 2017) was used to statistically estimate the relationship between the shear span-to-depth ratio and anisotropy ratio and the normalised strength by performing a polynomial regression analysis. The resulted regression formulae for the four-point bending test are presented in equation 8(8) $\mathit{\lambda }\left(\mathit{a}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.9099+0.0056{(\mathit{a}/d)}^2-0.0838\mathit{a}/\mathit{d}\pm 0.038E_L/E_T,\mathit{a}/\mathit{d}\le 8.2\\ 2.0465+0.0112{(\mathit{a}/d)}^2-0.2756\mathit{a}/\mathit{d}\pm 0.033E_L/E_T,8.2\le \mathit{a}/\mathit{d}\le 13.0\\ 2.3087+0.0069{(\mathit{a}/d)}^2-0.2453\mathit{a}/\mathit{d}\pm 0.021E_L/E_T,\mathit{a}/\mathit{d}\ge 13.0\end{array}\right.$(8) . Each formula consists of a second-order polynomial that describes a specific failure zone starting from lateral-torsional buckling, local buckling, and web crippling bounded by $\mathit{a}/\mathit{d}\ge 13.0$, $8.2\le \mathit{a}/\mathit{d}\le 13.0$, and $\mathit{a}/\mathit{d}\le 8.2$, respectively, when $\mathrm{L}/\mathit{d}=40.0$. The coefficient of determination R² that describes the accuracy of the regression model was 0.89. The anisotropy interaction with the a/d ratio was similar to its interaction with the L/d ratio. Increasing the $E_L/E_T$ ratio enhanced the lateral-torsional buckling resistance and decreased the size of web crippling and local buckling failure zones. In addition, the effect of laminate anisotropy was decreased when the a/d ratio increased to reach the lateral-torsional buckling zone.

Figure 9. Failure map showing the effect of shear span-to-depth ratio on I-shape pultruded FRP profiles subjected to four-point bending with different anisotropy ratios ($\mathit{L}/\mathit{d}=40.0$).

(8) $\mathit{\lambda }\left(\mathit{a}/\mathit{d},E_L/E_T\right)=\left\{\begin{array}{c}0.9099+0.0056{(\mathit{a}/d)}^2-0.0838\mathit{a}/\mathit{d}\pm 0.038E_L/E_T,\mathit{a}/\mathit{d}\le 8.2\\ 2.0465+0.0112{(\mathit{a}/d)}^2-0.2756\mathit{a}/\mathit{d}\pm 0.033E_L/E_T,8.2\le \mathit{a}/\mathit{d}\le 13.0\\ 2.3087+0.0069{(\mathit{a}/d)}^2-0.2453\mathit{a}/\mathit{d}\pm 0.021E_L/E_T,\mathit{a}/\mathit{d}\ge 13.0\end{array}\right.$(8)

It was noticed that the web crippling behaviour has approximately a constant function when the L/d ratio was changed (Figure 8) while it has a curved (second-order polynomial) function when the a/d ratio was changed. It was inferred that the web crippling strength is dominated by the a/d ratio that defines the distance between the supporting point and the nearest loading point. This distance controls the stress nature and distribution across the section. In web crippling, the stress of concern across the section is bearing stress caused by the transverse compression loading. Increasing the a/d ratio will shift this stress to bending stress with longitudinal compressive stress on the top flange and tensile stress on the bottom flange. Conversely, changing the L/d ratio, while maintaining a constant a/d ratio, will not change the stress distribution and will not have a significant effect on the web crippling behaviour. From Figures 8 and 9, it can be inferred that the lateral-torsional buckling and local buckling strengths are both significantly affected by the L/d and a/d ratios as evidenced by the steep $\mathit{\lambda }\left(\mathit{L}/\mathit{d}\right)$ or $\mathit{\lambda }\left(\mathit{a}/\mathit{d}\right)$ curves in their failure zones. Statistically, the percentile influence of the L/d and a/d ratios on the normalised strength $\mathit{\lambda }$ was estimated to be 61.8% and 38.2%, 53.5% and 46.5%, and 3.1% and 96.9% for the lateral-torsional buckling, local buckling, and web crippling failure modes, respectively.

Generally, the loading points move closer to the supporting points when the a/d ratio is decreased, which increases the beam’s tendency to fail by web crippling. On the contrary, when the a/d ratio is increased, the loading points move closer to the mid-span and the four-point bending becomes closer to the three-point bending. Consequently, the beam’s tendency to fail by local buckling or lateral-torsional buckling is increased depending on the L/d ratio.

The interaction between the L/d and a/d ratios was investigated to map the entire design space of the four-point bending test setup, as shown in Figure 10. This interaction was significant since the failure mode, as a function of the L/d ratio, was changed (shifted) depending on the a/d ratio. Equation 9(9) $0.0004(\mathit{L}/\mathit{d}{)}^2-0.0147\mathit{L}/\mathit{d}+0.6562>\mathit{\lambda }\left(\mathit{L}/\mathit{d}\right)\ge 0.0004(\mathit{L}/\mathit{d}{)}^2-0.0151\mathit{L}/\mathit{d}+0.4103\to \mathit{Localbuckling}$(9) describes the threshold boundaries between the different failure modes as resulted from the polynomial regression analysis (R² = 0.86). Each failure zone is bounded by a second-order polynomial that passes through all the inflection points (i.e. the points at which the slope suddenly changes) of all a/d curves. This design chart (Figure 10) can be used to predict the instability failure mode and its related buckling strength $\mathit{\lambda }$ based on the test setup parameters.

Figure 10. Failure map showing the interaction between the span-to-depth ratio and shear span-to-depth ratio of I-shape pultruded FRP profiles subjected to four-point bending.

$\mathit{\lambda }\left(\mathit{L}/\mathit{d}\right)<0.0004(\mathit{L}/\mathit{d}{)}^2-0.0151\mathit{L}/\mathit{d}+0.4103\to \mathit{Lateral}-\mathit{torsionalbuckling}$

(9) $0.0004(\mathit{L}/\mathit{d}{)}^2-0.0147\mathit{L}/\mathit{d}+0.6562>\mathit{\lambda }\left(\mathit{L}/\mathit{d}\right)\ge 0.0004(\mathit{L}/\mathit{d}{)}^2-0.0151\mathit{L}/\mathit{d}+0.4103\to \mathit{Localbuckling}$(9)

$\mathit{\lambda }\left(\mathit{L}/\mathit{d}\right)\ge 0.0004(\mathit{L}/\mathit{d}{)}^2-0.0147\mathit{L}/\mathit{d}+0.6562\to \mathit{Webcrippling}$

The failure mode shifts to lateral-torsional buckling ($\mathit{\lambda }\le 0.33$) when L/d ≥14.2 and $\mathit{a}/\mathit{d}\ge \frac{1}{4}\mathit{L}/\mathit{d}$. Whereas the failure mode transfers to local buckling ($0.33<\mathit{\lambda }\le 0.68$) when 36.5 > L/d ≥7.6 and $\mathit{a}/\mathit{d}\ge \frac{1}{6}\mathit{L}/\mathit{d}$. Increasing the a/d ratio reduces the L/d ratio required to attain lateral-torsional buckling or local buckling as the four-point bending transfers closer to three-point bending. The failure mode transfers to web crippling ($\mathit{\lambda }>0.68$) for any L/d ratio with $\mathit{a}/\mathit{d}\le \frac{1}{8}\mathit{L}/\mathit{d}$ or when L/d ≤26.4 and $\mathit{a}/\mathit{d}\le \frac{1}{6}\mathit{L}/\mathit{d}$ or when L/d ≤14.8 and $\mathit{a}/\mathit{d}\le \frac{1}{4}\mathit{L}/\mathit{d}$. Decreasing the a/d ratio increases the tendency to web crippling failure as it transfers the bending test setup closer to transverse bearing test setup subjected to transverse compression.

5. Conclusions

In this research, a numerical study was undertaken to investigate the structural instability failure modes of lateral-torsional buckling, local buckling, and web crippling of I-shape pultruded FRP profiles under bending. A novel modelling approach was established and verified against published experimental data and against the theory. The modelling approach introduces a Continuously Variable Step Size (CVSS) FE analysis to continuously examines and update the time step size based on artificially imposed variable geometric imperfection functions guided by scaled shape functions capturing the longitudinal and sectional deformations and instabilities. The validated modelling approach was used to perform parametric studies on the bending test setup parameters, namely the span-to-depth ratio L/d and shear span-to-depth ratio a/d. The effect of these parameters on the structural instability failure modes was investigated and quantified. The thresholds (i.e. the inflection points at which the failure mode changes to another one) between the instability failure modes were characterised graphically and statistically using failure maps and regression analysis. In addition, the interactions between the anisotropy ratio ($E_L/E_T$) of I-shape pultruded FRP profiles and test setup parameters were studied. These parametric studies aimed to introduce the appropriate test setup configurations to match the structural stability demands of real applications in construction. These parametric studies provide guidelines and recommendations for accurate testing and characterisation of the bending behaviour and structural instability failure of I-shape pultruded FRP beams. The following conclusions were summarised based on this study:

• The new modelling approach proved its capability and accuracy as a robust numerical tool that can capture different structural instability failure modes present under bending with low computational cost. The proposed modelling approach was faster by up to 33.0% compared to the default FEM analysis.

• When comparing three-point bending versus four-point bending, the lateral-torsional buckling, local buckling, and web crippling failure modes extended on 0.91, 0.88, and 1.14 of the L/d range under four-point bending compared to three-point bending, respectively. The related buckling strengths were higher under four-point bending by up to 12.6%, 7.2%, and 3.9%, respectively, compared to three-point bending for the same L/d ratio.

• The web crippling strength is dominated by the a/d ratio while the lateral-torsional buckling and local buckling strengths are both significantly affected by the L/d and a/d ratios. The percentile influence of the L/d and a/d ratios on the normalised buckling strength ($\mathit{\lambda }$) was estimated to be 61.8% and 38.2%, 53.5% and 46.5%, and 3.1% and 96.9% for the lateral-torsional buckling, local buckling, and web crippling failure modes, respectively.

• The laminate anisotropy ratio $E_L/E_T$ contributed up to 11.9%, 35.3%, and 46.2% of the lateral-torsional buckling, local buckling, and web crippling strengths. The effect of laminate anisotropy decreased when the L/d ratio increased to reach the lateral-torsional buckling zone as the beam span length has a higher effect compared to the laminate layup properties.

• The interactions between the L/d and a/d ratios in four-point bending were significant in shifting the failure modes. The failure mode shifts to lateral-torsional buckling ($\mathit{\lambda }\le 0.33$) when L/d ≥14.2 and $\mathit{a}/\mathit{d}\ge \frac{1}{4}\mathit{L}/\mathit{d}$. Whereas the failure mode transfers to local buckling ($0.33<\mathit{\lambda }\le 0.68$) when 36.5 > L/d ≥7.6 and $\mathit{a}/\mathit{d}\ge \frac{1}{6}\mathit{L}/\mathit{d}$. Increasing the a/d ratio reduces the L/d ratio required to attain lateral-torsional buckling or local buckling. The failure mode transfers to web crippling ($\mathit{\lambda }>0.68$) for any L/d ratio with $\mathit{a}/\mathit{d}\le \frac{1}{8}\mathit{L}/\mathit{d}$, or when L/d ≤26.4 and $\mathit{a}/\mathit{d}\le \frac{1}{6}\mathit{L}/\mathit{d}$, or when L/d ≤14.8 and $\mathit{a}/\mathit{d}\le \frac{1}{4}\mathit{L}/\mathit{d}$. It is recommended to consider these guidelines when designing test setups for accurate characterisation of structural instability failure modes.

The findings in this research are limited to I-shape cross-sectional geometry subjected to bending loading. Future studies that consider other cross-sectional shapes and loading conditions are needed in order to develop the current modelling approach to include other structural instability failure modes. Such studies will enhance the understanding and provide comprehensive recommendations and design guidelines for different types of test setups and failure modes.

Acknowledgments

The author gratefully acknowledges the technical support provided by Tafila Technical University to complete this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data presented in this study are available on request from the corresponding author.

Additional information

Funding

No funding was received to assist with the preparation of this work.

Notes on contributors

Mohammad Alhawamdeh

Mohammad Alhawamdeh is an assistant professor at Tafila Technical University in Jordan. He is a specialised researcher in numerical simulation and optimisation of composites. Alhawamdeh holds Ph.D. (2022) and M.Sc. (2016) in structural engineering. He is a skilled designer using computer-aided software and programming languages, such as Abaqus, Python, and MATLAB. Alhawamdeh is a motivated researcher who interacts with analytical, experimental, and numerical approaches. His field of research includes composites analysis and design, numerical optimisation of composites, and finite element modelling. His work aims to develop advanced composite pultrusion for civil infrastructure applications and enable rapid commercialisation of pultruded profiles in this sector.

Appendix A.

List of notations

ѱ: geometric imperfection function of the profile (mm)

ɸ_{b, f, w}: shape function of the beam, flange, or web, respectively (mm)

ɣ_{b, f, w}: shape function weight of the beam, flange, or web, respectively (%)

δ_{b, f, w}: normalised displacement of the beam, flange, or web, respectively

$\mathit{\lambda }$: normalised buckling strength

$v_{\mathit{LT}}$: in-plane Poisson’s ratio of the profile

a/d: shear span-to-depth ratio of the beam

d: cross-sectional height of the beam (mm)

E_L: longitudinal elastic modulus of the profile (MPa)

E_T: transverse elastic modulus of the profile (MPa)

G_LT: in-plane shear modulus of the profile (MPa)

L/d: span-to-depth ratio of the beam

L: span length of the beam (mm)