The first experimental research stage focuses on fiber quantification in the GFRP composites. The standard ISO 1172:1996 “Textile-glass-reinforced plastics—Prepregs, molding compounds and laminates—Determination of the textile-glass and mineral-filler content—Calcination methods” defines the test procedure requirements. An electric furnace Snol 7.2/1100 (Snol, Lithuania) was used for the heating tests. For that purpose, small fragments (≈10 g), cut from the webs and corners of the profiles (Figure 2), were dried until the constant weight at 105 ± 5°C (the drying time = 24 h). After that, the specimens were heated at 625 ± 5°C until the polymer evaporated. The heating regime: 2 h to reach a temperature of 625°C and 3 h to achieve a constant weight. The heating procedure determines the mass percentage of fibers w _f , and the following equation defines the volume fraction V _f : V f = w f ⋅ ρ m w f ⋅ ρ m + ( 1 − w f ) ⋅ ρ f , (1) where ρ _f and ρ _m are the fiber and matrix densities (Table 1). The above parameter is necessary for the numerical modeling to describe the fiber reinforcement area in the cross-section (Gribniak et al., 2012, 2021). However, this testing methodology gets only an approximate estimate of w _f because of neglecting the heat-resistant aggregates’ effect (Bazli and Abolfazli, 2020).

Figure 3 shows the analysis results of three (#1–#3) web and corner samples for each profile type. This figure demonstrates that the estimated fiber content (71.0%) exceeds the mass fraction declared by Fiberline (60 ± 5%). The neglecting of the weight of the heat-resistant aggregates cannot explain such a big difference that can be a consequence of the complex reinforcement architecture if the manufacturer specifies only longitudinal fiber content. Fiberline (2003) states that various types of roving and intricate weaves and mats compose the reinforcement system of the pultruded structural profiles. Unfortunately, the detailed information about the fiber percentage and orientation is missed in the Fiberline datasheets.

On the contrary, the average Ugira fiber fraction (75.9%) agrees with the declared value (75%), indirectly indicating the simplicity of the reinforcement layout. Eq. 1 results in the average volume fraction of 56.2 and 67.8% for Fiberline and Ugira profiles. The calculations employ density values specified in Table 1.

Figure 3 also demonstrates almost identical fiber content in the web and corner specimens for both profiles. It means that the current pultrusion technologies solved the resin localization problems at the web-flange junction identified previously (Feo et al., 2013).

The tensile coupon tests were carried out for Ugira profiles within the framework of this experimental program by following the ASTM D 3039/D 3039M-17 “Standard Test Method for Tensile Properties of Polymer Matrix Composite Materials” requirements. The mechanical properties, i.e., the tensile strength and elastic modulus, were determined in the pultrusion direction, corresponding to the filament orientation. All coupons were tested using a 250 kN capacity servo-hydraulic machine Z250 (Zwick, Germany). Four coupons were tested for Ugira profiles. Table 2 summarizes the tensile test results. This table also includes the tensile test outcomes provided by Fiberline (2021).

The tensile coupon tests for the Fiberline profiles were not carried out in this illustrative study. Thus, this manuscript provides a general view of the problems related to developing a reliable numerical model regarding the incomplete information provided in the producers’ certificates.

The microstructure analysis is a part of the current experimental study. The filaments after the heating tests (see section “Fiber Content Analysis” above) and the web and corner microsection specimens of the SHS profiles (Figure 2) define the investigation objects. A field emission scanning electron microscope JSM-7600F (Jeol, Japan) was used for the microstructure analysis. The scanning electron microscopy (SEM) parameters were as follows: 10 kV voltage, 11–13 mm distance to the specimen surface and magnification varied from ×50 to ×1200.

Figure 4 shows the SEM images of separate fibers obtained during the heating tests. The difference in fiber diameter is the center outcome of these images—the Fiberline samples vary in the 16–22 μm range; the alternative fibers’ diameter varies from 32 to 38 μm.

Abrasive grinding and thin Au layer electrically conductive covering in a vacuum with a Quorum Q150R ES equipment (Quorum, United Kingdom) prepared the microsection specimens. Figures 5, 6 show the SEM images of the web and corner specimens’ surfaces normal to the pultrusion direction. The Fiberline samples demonstrate a more dense structure than the Ugira specimens do. Independently of the sampling (Figures 5A,B, 6A,B), the fiber and matrix contact zones have no visible defects. On the contrary, deficiencies and voids of size varying from 10 to 33 μm are characteristic of the Ugira samples (Figures 5C,D, 6C,D).

The material heterogeneity and anisotropy of GFRP composites alter material performance within the pultruded cross-section (Correia, 2013). Besides, a complex stress-strain behavior contributes to the tensile and compressive material properties distinction in a flexural sample (Insausti et al., 2020). Thus, three-point bending tests were included in this experimental program, determining the mechanical properties of structural profiles. This study employs three Fiberline (50 × 50 × 5 mm) and five Ugira (40 × 40 × 3.2 mm) SHS profile samples. All the test samples had a length of 1000 mm and were tested until failure. Steel rollers and bearing steel plates protected the GFRP profile from the local damage. Figure 7 shows the loading setup.

A 5 MN capacity servo-hydraulic machine LFV 5000 (w+b, Switzerland) loaded the profiles under the displacement control, measuring the applied load with a 100 kN capacity load cell. A digital image correlation (DIC) system captured surface deformation of the profiles; paired 100 mm range and 0.01 mm precision linear variable displacement transducers (LVDT) measured the mid-span vertical displacements. An Almemo 2890-9 data logger acquired the test data every second.

Figure 8 shows the average mid-span displacements measured by the paired LVDTs (Figure 7A). The load-displacement diagrams of all specimens (i.e., S _Fl 1-S _Fl 3 and S _Ug 1-S _Ug 5) are almost linear until the ultimate load. The samples also revealed similar flexural stiffness and failure mechanism—crushing of the compressive zone followed by cross-section integrity loss. Figure 9 shows the typical failure cases.

The DIC system (Figure 7B) captures the relative displacements of the surface points (pixels) recognized at the monitoring surface at the rate of one image per 500 N with the aid of the Davis 8.1.6 software by Lavision. That allows tracking relative movements of any pixels chosen at the image processing stage (Misiunaite et al., 2020).

The FE model employs the smeared reinforcement concept using commercial FE software Atena (Cervenka, 2002). Figure 10 sketches the FE modeling concept following that virtual filaments are uniformly distributed inside solid finite elements. The straight fibers have no specific geometry—the volume percentage and orientation (regarding global coordinate system) describe the mechanical properties of the reinforcement having a perfect contact to the matrix (i.e., no slip between the filament and finite elements is allowed). The fibers are not resistant to axial compression. The solid finite elements can include several reinforcement systems (layers) smeared in different directions to represent a complex filaments architecture characteristic of particular FRP pultruded components.

In this study, the 3D solid finite elements describe the polymer matrix, assuming the fracture mechanic principles for the tensile failure and the plasticity approach to the compressive failure. Figure 11 shows the softening law in compression described with a linear descending branch of the diagram. In this figure, f _c is the compressive strength; f _c0 is the onset of non-linear compressive behavior; ε _cp is the plastic strain at the compressive strength; ε _d is the compressive ductility. The ratio between the plastic displacement and crack band size defines the strain ε _d . The crack band determines a fictitious crushing zone that is related to the finite element size. The principal compressive stresses describe the failure plane normal to the stress direction; the post-peak compressive strain localizes in this plane. Such formulation reduces the FE mesh dependency of the model (Cervenka, 2002).

A layer of the smeared reinforcement models the glass filaments. In such a way, a separate 1D material model, coinciding with the pultrusion direction, describes the reinforcement. A perfectly elastic-brittle constitutive law defines the tension failure of the fibers, which do not resist the compression stresses.

Figure 12 shows the numerical model built using 3D isoparametric brick eight-node finite elements with eight integration points and an average 3 mm size. The modeled steel plates protect the GFRP profile from stress concentration at the supports and load application point. These boundary and loading conditions (Figure 12A) correspond to the physical tests (Figure 7A). The FE model of the Ugira profile (Figure 12B) reflects the shape imperfection (i.e., a relative rotation of the inner rectangular shape in the cross-section, Figure 12C) that resulted from the manufacturing flaw. Several FE simulations were carried out, varying only the assumed fiber content Vf. Table 3 summarizes assumed parameters of the numerical models. The following section describes the Vf quantification principles.

Tables 1, 2 provide experimental values of all material parameters necessary for the numerical modeling of the Ugira profile. Therefore, simulations of this profile illustrate the proposed fiber efficiency concept and the limitations of the standard testing procedures. The material parameters listed in Table 1 describe the elastic constitutive models of the vinylester resin and smeared reinforcement in all simulations. The equal compressive and tensile strength values were assumed for the polymer.

Model 1 accepts the volume percentage V _f = 67.8%, corresponding to the fraction experimentally determined in the “Fiber Content Analysis” section above. The tensile coupon test (section “Mechanical Properties of GFRP Composite” above) defines the fiber content of Model 2—the external tensile load P induces two internal forces acting on the polymer matrix N _m and glass fibers N _f : P = N m + N f . (2)

The average axial strains can express the above force components as follows E c ⋅ A c ⋅ ε c = ∑ i = { m , f } E i ⋅ A i ⋅ ε i , (3) where E, A, and ε are the elastic modulus, area, and average axial strain; the indexes c, m, and f correspond to the GFRP composite, polymer matrix, and glass filaments. Note, the volume content is the center parameter describing the reinforcement area A _f ; the fiber mass percentage is not suitable for this purpose.

The Navier-Bernoulli hypothesis (planar strain distribution within a cross-section height) postulates equality of the average strains: ε c = ε m = ε f , (4) simplifying Eq. 3. Thus, the following formula expresses the relationship between the volumetric fiber fraction V _f and the areas A presented in Eq. 3: A m = A c ⋅ ( 1 − V f ) ; A f = A c ⋅ V f . (5)

The solution of Eq. 3, assuming the simplifications of Eqs 4, 5 and the experimental value of the composite coupon elastic modulus E _c,exp , determines the efficient reinforcement fraction V _f,eff as follows V f , e f f = E c , e x p − E m E f − E m . (6)

Tables 1, 2 describe the elastic moduli for determining the efficient volume V _f,eff = 53.2% assumed in Model 2. This calculated fraction is well below the experimentally determined value of 67.8% (taken in Model 1).

The vertical displacement, applied in small increments (0.3 mm), determined the load-bearing capacity of the models. The deformation problem was solved in a 3D formulation by the Newton-Raphson iteration procedure. Figure 13 demonstrates the inadequacy of the volume fraction V _f,eff determined from the tensile coupon test results to predict the deformation behavior of the GFRP profile (Model 2). Besides, the fiber fraction was tailored to represent adequately the experimental diagrams shown in Figure 8B. The iterative trial-and-error process resulted in the content V _f,eff = 62.0% assumed in Model 3; Figure 13 demonstrates the deformation prediction adequacy.

Misiunaite et al. (2020) described the verification procedure principles. This manuscript compares the local deformations predicted by Model 3 and the DIC surface monitoring results. The DIC system (Figure 7B) estimates local strains near the steel plate edge, using the 10 mm virtual strain gauges, as Figure 14A shows. Model 3 determines the relative displacements of the nodes corresponding to the virtual gauge positions (Figure 14B).

Figure 15A shows the strain distribution over the profile height identified by the DIC system capable of capturing relative pixel displacements belonging to the vertical monitoring surface (colored areas in Figure 15). Thus, the DIC diagrams (Figure 15A) do not reach the profile horizontal surfaces. The rounded corners of the SHS profile (Figure 12C) define the unidentified areas. The FE model did not represent the corners’ roundness (Figure 12B). That is a consequence of the engineering simplifications, determining the FE modeling principles. Therefore, the simulated diagrams (Figure 15B) have no limitations characteristic of the DIC system. The strain profiles are related to the corresponding external load P.

The DIC and numerical simulation results’ comparison (Figure 15B) proves the above inference—Model 3 correctly represents the compressive strain increase while the smeared reinforcement does not resist compression (Cervenka, 2002). The comparative analysis of these diagrams proves sufficient adequacy of Model 3 predicting non-linear deformation behavior of the GFRP composite profile. In most cases, the differences between the tensile strains calculated using FE and estimated by the DIC system do not exceed 10%. The average strain differences do not exceed 15%, with the maximum errors located near the neutral axes (where the strain magnitudes are comparable with the measurement tolerance). Failure of the compressive zone governs the predicted ultimate load (Figure 13). That agrees with the results shown in Figure 9 and the literature findings (Correia et al., 2011; Wu and Bai, 2014; Wu et al., 2019; Alaedini et al., 2021; Almeida-Fernandes et al., 2021).

The engineering application defines the optimum FE mesh requirement—a sufficiently coarse mesh must ensure acceptable prediction accuracy. Using the 1.5 and 1.0 mm average FE size, two additional simulations verified the Model 3 prediction adequacy. (This manuscript focusing on the physical investigation aspects does not present those simulation results.) The mesh sensitivity analysis revealed the following outcomes.

1) The mesh size does not affect the flexural stiffness prediction results. That is the expected consequence of the smeared reinforcement model’s robustness for engineering application (Cervenka, 2002; Gribniak et al., 2010).

The absence of materials parameters (Table 1) essential for the proposed FE approach makes the Fiberline profile analysis only illustrative. Numerical simulations of the bending profile employed the same procedure as specified in the “Finite Element Model” section; the “Flexural Tests” section describes the modeling object. Figure 16 shows the corresponding numerical model built using 3D isoparametric brick eight-node finite elements with eight integration points and an average size of 5 mm. The modeled steel plates preserve the GFRP SHS profile from stress concentration at the supports and load application point. Figure 16A defines the boundary and loading conditions, representing the physical tests (Figure 7A). Figures 16B,C show that the Fiberline model has no geometry imperfections characteristic of the Ugira profile (Figure 12B). However, both FE models ignore the corner roundness of the SHS profiles (Figures 12B, 16B) because of the engineering simplification.

In the same manner, as described in the “Finite Element Model” section, the 3D solid finite elements describe the polymer matrix; smeared reinforcement represents the glass filaments (Figure 10). The respective constitutive model of the GFRP material employs the polyester resin and E-glass fiber properties obtained from the literature (El-Wazery et al., 2017; Fiberline, 2021) and presented in Table 1.

Two FE simulations were carried out, assuming different fiber volume fractions. Model 4 considers a profile with the 56.2% filament content estimated from the mass fraction obtained in the “Fiber Content Analysis” section. Model 5 assumes the fiber volume fraction of 44.0% that corresponds to the manufacturer specified nominal mass fraction (60%) transformed by Eq. 1 and the constituent material densities ρ _m and ρ _f defined by Fiberline (2021).

Figure 17 compares the numerical prediction and the experimentally measured vertical displacements at the mid-span. The FE simulations consider the load increments applied in small (0.4 mm) steps. Model 4 demonstrates the inadequacy of the filament content estimated using the mass fraction values obtained during the heating tests (the “Fiber Content Analysis” section). Still, Model 5 demonstrates good agreement with experimental observations.

The standard tests consider material fragments. That is acceptable for homogenous media such as metals and pure polymers. However, the fiber reinforcement makes the coupons cut from FRP composites unrepresentative for investigating material properties. Figure 18 shows a microstructure corresponding to the tension coupon cut-side surface consistent with the pultrusion direction of the Ugira profile (the “Mechanical Properties of GFRP Composite” section). Numerous damaged filaments result from such a cut. This figure also demonstrates the bond damages of the separate fibers (Figure 18B). Similar surfaces define the boundaries of the coupon samples for determining the effective fiber fraction V _f,eff in Model 2 (the “Fiber Efficiency Analysis” section). The filament damages (Figure 18A) explain the reinforcing efficiency losses (compare V _f,eff of Model 2 and Model 3).

The numerical analysis (the “Analysis of the Fiberline Profile” section) of the Fiberline composite elastic modulus (Table 2) highlights a remarkable underestimation of the fiber reinforcement efficiency based on the tension coupon tests (Fiberline, 2021). Under the assumption of the elastic moduli of the composite and constituents from Tables 1, 2, Eq. 6 defines V _f,eff = 26.9% that is well below the fiber fraction 44.0% assumed in Model-5 (Figure 17). That indicates the inapplicability of the standard coupon tests for constitutive modeling of FRP composites. Furthermore, the material properties estimation error increases with the decrease of the sample size, i.e., with the cut-side area increase regarding the sample cross-section dimensions.

This study also questions the clarity of the manufacturer datasheets. The complex internal architecture of the Fiberline composite is a possible explanation of the inconsistency between the declared (≈60%) and assessed in the “Fiber Content Analysis” section (71.0%) fiber contents. Figures 6A,B illustrate the differences in the filament orientation regarding the pultrusion direction, but the producer does not provide any detailed information on the filament fraction quantification. The FRP constituencies (fibers and polymer) characterization results can help identifying the effective fiber fraction V _f,eff essential for the numerical simulations.

Table 4 describes the numerical predictions’ accuracy described in the flexural stiffness EI and ultimate load P _u terms. The first three models represent the Ugira profile, and the Fiberline SHS sample was the modeling object of the remaining two simulations. The identical profile models had the same geometry, FE mesh, and material parameters; the volumetric fiber content was only the difference. The elastic analysis defines the flexural stiffness as follows E I = P ⋅ l 3 48 ⋅ u , (7) where P and u define the load and corresponding vertical displacement; l is the loading span (= 0.9 m, Figures 7A, 16A). The elastic behavior of the flexural profiles (Figure 8) makes the above simplification possible. The stiffness EI was calculated for the experimentally determined P-u pairs at all ascending load stages. Table 4 shows the EI and P _u values averaged for the identical profiles. The difference between the predicted and experimental outcomes divided by the test value defines the prediction error.

The modeling results (Table 4) demonstrate that the assumed fiber volume V _f affects the predictions. Accounting for the variation of the experimental outcomes, the significant errors (in the statistical sense) are characteristic of the Ugira profile stiffness results. Furthermore, the prediction error exceeds the V _f assessment error (compare the EI results of Model 1 and Model 2 with the reference Model 3). The analysis of the Fiberline predictions identifies a similar trend—inadequate V _f assumption increases the EI prediction error from 3.5% (Model 5) to 21.8% (Model 4). However, the opposite tendency in the ultimate load prediction exists—the fiber content increase remedies the P _u predictions. That can be a consequence of the effect of the additional fibers identified during the heating tests (the “Fiber Content Analysis” section) regarding the manufacturer specified content.

An efficient balance between the prediction accuracy and computation time characterizes the developed FE approach that does not require specific descriptions of reinforcement geometry and refined meshes necessary for modeling the discrete fibers (Garnevičius et al., 2020). The Atena software also allows smearing the reinforcement in different directions to represent the complex fiber architecture (the “Finite Element Model” section). Unfortunately, the manufacturer provides no sufficient information to develop such a model. Therefore, the updated datasheets must include information about the fiber percentage and orientation in different directions. The proposed modeling methodology can also determine the fiber efficiency in reinforcing the polymer matrix useful for developing efficient structural components.

The pultrusion direction governs the mechanical resistance of FRP profiles (Correia, 2013). In three-point bending, the center load stresses the profile in both pultrusion and transverse directions. Under such loading conditions, the analogs metallic SHS profiles face the web crippling because of the web yielding or web buckling failure mechanisms. The failure mode depends on the relative web stiffness described in terms of the slenderness factor (Misiunaite et al., 2020): λ = b t f y / E , (8) where b and t are the width and thickness of the flat web part; f _y and E are the yield strength and elastic modulus of the material.

For isotropic materials, the principal stress and principal strain orientations coincide at a 45° angle. This condition is not valid for orthotropic FRP profiles, where the pultrusion direction does not describe the principle stress-strain distribution (Insausti et al., 2020). Moreover, a low shear modulus of FRP composites (regarding the steel) causes intensive shear deformations of short beams also known, as shearing strain (Peng et al., 2020). The excessive shear strains define a central problem in pultruded FRP beams because of the minimal contribution of longitudinal fibers (Correia et al., 2011; Wu et al., 2019; Almeida-Fernandes et al., 2021). The misbalance in longitudinal and transverse material properties leads to non-proportional angular distortion and shearing failure of the web, shown in Figure 9. Thus, the web crippling mechanism of the pultruded GFRP SHS profile combines the web-flange junction’s separation and the web shearing/crushing.

The strain distribution within the cross-section height of the GFRP SHS beam is non-linear with the increase of inelastic deformations. Moreover, the Navier-Bernoulli hypothesis is not valid at already low loading levels (e.g., 3 kN, Figure 15A). That is the consequence of the synergetic effect of the tensile and compressive elastic modulus distinction in the GFRP composite (Insausti et al., 2020), and the slender filaments instability in compression causes the onset of buckling in the compression flange (Correia et al., 2011; Misiunaite et al., 2020). However, the proposed approach to FE modeling based on the smeared reinforcement concept adequately covers both above material aspects characteristic of the pultruded GFRP composites.

The Fiberline fibers have uneven rough surfaces with small particles having an approximate 1–3 μm size (Figures 4A,B) that can potentially improve the bond performance with a polymer matrix (Papon et al., 2020). On the contrary, an untreated glassily surface is characteristic of the Ugira samples (Figures 4C,D). That defines the object for further modifications.

Figures 5, 6 show the SEM images of the web and corner microsections normal to the pultrusion direction. The reference Fiberline specimens demonstrate a more dense microstructure without visible defects of the fiber contact zones with polymer matrix than the Ugira samples do. With the increased fiber diameter, the interactional defects observed in the alternative specimens (Figures 5C,D, 6C,D) can further affect the bond performance and reduce the efficient fiber fraction V _f,eff of the Ugira profiles identified in the “Fiber Efficiency Analysis” section. The bonding mechanisms identified in reinforced concrete composites (Gribniak et al., 2016, 2020; Rimkus and Gribniak, 2017) support the above inference—the same reinforcement area in a higher number of bars (with reduced bar diameter) can substantially improve the bond performance of the composite material. Thus, the fiber diameter reduction for the same fiber-reinforcement fraction can potentially increase its efficiency.