3 Proposed theoretical approximation Equation Chapter Se ction
3.2 theoretical approximation
3.2.2 Numerical approximation
The influence of the fillets could be added to extend the theoretic approximation. However, the actual influence of fillets was difficult to describe with simple mechanics, as the vector of xz changes direction while the flux is maintained. In this chapter two methods are proposed, in which the vector is remained in plane with the web. Method 1 is based on the principle of a constant shear flow in the web; the shear stress is reduced directly proportional to the increasing width at the fillets. Method 2 uses Jourawski’s formula (3.1) to obtain a shear stress distribution for a discontinuous width.
Incorporating the fillets introduces a discontinuous function for the width. Therefore, both methods will consist of summations for different sections. To apply Jourawski’s formula, the first moment of inertia at the fillets needs to be determined. The first moment of inertia may be defined through equation (3.9).
2
2 2 1
2( ) 3 ( ) cos 2
fillet 4 3
z r r z
S z z r r r h rh h
r z r
(3.9)
Jourawski’s formula is known to be less accurate for non-rectangular shapes (Carpinteri, 1997).
Moreover, a constant shear flow is preferred, therefore method 1 is preferred. Considering approach 1, a decrease in shear stress due to the increasing width allows an increase in axial stress using the Von Mises equation. Formulas are obtained for the development of shear- and axial stresses in the range V/Vpl=0 to 1. These formulas are subdivided to their corresponding sections of the web, as the width is discontinuous. Subsequently, these stresses are integrated over the height of their sections for increasing shear stresses, and create sets of MV values. These sets can be curve fitted to obtain M-V interaction formulas, however, they will be section dependent. The datasets of HEA, HEAA, HEB, HEM and IPE are created with loops in Maple and added in APPENDIX VII, a worksheet for an individual section is added in APPENDIX VI. The datasets indicate HE sections are less affected by shear stress than IPE sections, since their moment capacities are less dependent on the web. Moreover, the IPE datasets, as illustrated in Figure 15, seem to indicate that the influence of shear stresses on the moment capacity increases with height. Therefore, M-V interaction design rules might need to incorporate such variables, still minding the usability of the design rule.
4 N UMERICAL MODEL
EQUATION CHAPTER 4 SECTIO N 4To ensure the numerical model accurately describes the mechanical behavior several modelling aspects were evaluated in a parametric study. The influence of the modelling aspects on the behavior of the model were evaluated. Subsequently, to validate the constitutive behavior of the Abaqus model the displacement and rotations were compared to corresponding functions that were derived from mechanics. To extend the validation in plasticity and instability, the numerical model was compared to experimental tests.
4.1
P
ARAMETRIC STUDY 4.1.1 Geometry and setupConsidering the purpose of this research the IPE section was chosen as the basic type section, as it is most prone to instability. During validation, no deviation was made from existing sections. Two lengths were selected to test all aspects for shear dominated(SD)as well as bending dominated(BD) cases. The shear dominated case dimensions were based on the threshold for Euler Bernoulli beam theory.
Eurocode 2-5.3.1(3) does not distinguish aspect ratios for deep beams in 3 or 4-point bending test, instead, it merely states a requirement of 3h between the supports. Therefore, the definition (a=2h) for 4-point bending tests as is illustrated in Figure 16 was obtained from ACI 318 section 10.7. The bending dominated case was defined as (a=9h). During validation, a true stress strain relationship including strain hardening was used as is illustrated in APPENDIX VIII.I.
NEN-EN 1993-1-5.1(2) states webs with hw/t < 72ε should have transverse stiffeners at the supports, since this is well below the slenderness criteria for class 3 illustrated in Figure 2, reference nodes 3/6 and 5/8 were used to model rigid bodies as stiffeners.
The middle lines of the flanges were constrained in X direction, to prevent lateral torsional buckling.
Sections S2 and S3 were used to obtain the bending moment during loading. RP10 was used to read the displacement at mid-span. RP1, RP2 were used to apply a bending moment to the beam, they were coupled to sections S1, S4 respectively through kinematic coupling. RP10 and RP11 were used to apply loads by means of an equation coupling, as will be explained in depth in chapter 4.1.7.
Figure 16-Geometry and setup numerical model parametric study
Furthermore, additional abbreviations were used to address variables in graphs and tables. For instance, the number of elements(EL), force controlled (FC) or displacement controlled(DC), or if half (HALF)a model with symmetry axis was used, otherwise they were specified in the chapter itself.
4.1.2 Elements
According to Dekker, (Dekker, 2018) the fillets influence the shear stress distribution, and
consequently, the shear resisting capacity of steel sections. Therefore, the fillets were not modelled as a simplification of reality and had to be modeled with solid elements. Linear brick elements, except from the C3D8I element, are prone to hourglassing and shearlocking. Therefore, the linear brick elements were compared to the quadratic brick element that should behave well in this situation. The objective was to find an element that behaves well with the least computational power. An overview of the FEM-analyses in this chapter is presented in Table 1, results of the total test population are presented in APPENDIX IX.II.
Table 1- Overview of simulations for evaluating element types
The graphs in APPENDIX IX.II indicate that shear dominated cases were not affected by the type of element, however, bending dominated cases were. Both linear elements C3D8 and C3D8R appeared to be overly stiff in plasticity. Overly stiff behavior was to be expected from the C3D8 element due to shearlocking, the additional capacity while using the C3D8R element in Figure 17 was explained differently. The C3D8R element showed excessive displacements, probably due to hourglassing.
Consequently, large strains correspond to high stresses due to strain hardening. The C3D8R element widely overshoots the plastic bending moment capacity, and thereby confirms strain hardening occurred. Therefore, both C3D8R and C3D8 elements were not suitable elements for the numeric model. The C3D8I and C3D20R element behaved similarly, and accurately obtained the plastic bending moment capacity. Therefore, the C3D8I element was chosen as it required less computational power.
Figure 17- M-δ graphs of 4-point bending test with linear and quadratic elements 0
0 200 400 600 800 1000 1200
M/Mpl[-]
Bending dominated Shear Dominated
C3D8 X X
C3D8R X X
C3D8I X X
C3D20 X
C3D20R X X
4.1.3 Symmetry axis
Addition of a symmetry axes could reduce the size of the model in half, decreasing the required computational power. The symmetry axis was placed at mid span, by constraining translations in Z and X directions. An overview of the FEM-analyses in this chapter is presented in Table 12, the results of the total test population are presented in APPENDIX IX.III.
Table 2- Overview of simulations: Symmetry axis
Without symmetry axis With symmetry axis
FC DC FC DC
4 elements (EL) X X X X
2 elements (EL) X X X X
The buckling mode without symmetry axis indicates the compression flange is angled at mid-span.
However, the symmetry axis pins nodes horizontally and disables rotations as is showed in Figure 18.
(a) Eigenvalue: 376,94 (b) Eigenvalue: 399,8
The symmetry axis affected the buckling modes and increased the eigenvalues, making it an unsuitable method to model a 4-point bending test. The graphs in APPENDIX IX.III support this conclusion. The python input file created an image for every increment that was related to a maximum in the force-displacement(F-δ) graph. The F-δ graphs without symmetry axis as illustrated in Figure 19 , showed two maxima, in which the images clearly showed a first, and second sine wave in the compression flange. These maxima were not present in the F-δ graphs with symmetry axis.
Moreover, the obtained bending moment capacity of the model with symmetry axis was higher, due to the higher eigenmodes. Therefore, a symmetry axis was not applicable in this configuration.
0 10000 20000 30000 40000 50000 60000
0 200 400 600 800 1000 1200
Force[N]
Displacement[mm]
4EL-DC-BD 4EL-FC-BD 4EL-HALF-DC 4EL-HALF-FC Figure 18- Instability with (b) and without (a) symmetry axis
4.1.4 Mesh size section
Prior studies (Dekker, 2018) indicated that shear stresses distributed outwards over the inner sides of the flanges, therefore, a sufficient number of elements had to be used over the thickness of the flanges and web. The number of elements was decreased from 8 to 4, and subsequently 2 as is illustrated in Figure 20. An overview of the FEM-analyses is presented in Table 3, the results of the total test population are presented in APPENDIX IX.IV.
Table 3- Overview simulations: Number of elements over thickness
The computation time of 8 elements BD was excessive, and therefore, this configuration was modeled with a symmetry axis. As explained previously, the validity of the model with symmetry axis expires after reaching Mpl. The graphs in APPENDIX IX.IV, and Figure 21 indicate no conclusive effect of the amount of elements on the results is apparent. The Abaqus manual recommends a minimum of 4 elements over the thickness of plates, therefore, this minimum was adopted for further analyses.
0 10000 20000 30000 40000 50000 60000
0 200 400 600 800 1000 1200
Force[N]
Displacement[mm]
8EL-HALF-DC-BD 2EL-DC-BD 4EL-DC-BD
Bending dominated Shear dominated
FC DC FC DC
2 EL X X X X
4 EL X X X X
8 EL X (HALF) X(HALF) X X
tf: 2 elements tf: 4 elements tf: 8 elements tw: 4 elements tw: 4 elements tw: 8 elements
Figure 20- Mesh study on the amount of elements over thickness of plates
4.1.5 Mesh size aspect ratio
Over the length of shear- axial stress interaction (section a) the aspect ratio was varied. In literature, many different aspect ratio limits are prescribed, however, the correct aspect ratio of an element is very context dependent. Since instability and large deformations were expected, high aspect ratios should be avoided. A well-regarded safe aspect ratio of 1:10 is used in the regions near loads and supports.
Figure 22- regions of the beam with varying mesh size aspect ratios
To avoid erratic behavior the transition of aspect ratios between adjacent regions was smoothened.
The Abaqus manual 6.10 prescribed an aspect ratio limit of 1:20, and therefore, this ratio constituted the limit of this aspect ratio study. Subsequently, the aspect ratio was decreased to 1:10, and 1:5. An overview of the FEM-analyses is presented in Table 4, the results of the total test population are presented in APPENDIX IX.V.
Table 4- Overview simulations: Longitudinal mesh size aspect ratio
Mesh size aspect ratio Bending dominated Shear dominated
1:5 X X
1:10 X X
1:20 X X
The graphs in APPENDIX IX.V, and Figure 23 seem to indicate no conclusive effect of the aspect ratio on the behavior is apparent. Even though the aspect ratio 1:20 seemed to behave well in this configuration, the aspect ratio of 1:10 was adopted to avoid potential problems, while a reasonable computation time was maintained.
0 10000 20000 30000 40000 50000 60000
0 200 400 600 800 1000 1200
Force[N]
Displacement[mm]
1:5-LONG-BD 1:20-LONG-BD 1:10-LONG-BD
4.1.6 Force and displacement controlled Loading
The beam may be loaded by either an imposed displacement or force. An overview of the analyses is presented in Table 5, the results of the total test population are presented in APPENDIX IX.I.
Table 5- Simulation overview: Force controlled or. Displacement controlled
Force controlled or displacement controlled loading should not differ since the RIKS method solves both load and displacement simultaneously by computing the arc length. In contrary to the static general step, forces and displacement are not required to keep increasing while using Riks.
However, as can be seen from the graphs in APPENDIX IX.I, and Figure 24, displacement controlled analyses seemed to be more stable as displacements keep increasing and forces do not. Therefore, displacement controlled analyses were preferred.
Figure 24- F-δ graph of 4-point bending test loaded displacement controlled or force controlled 0
10000 20000 30000 40000 50000 60000
0 200 400 600 800 1000 1200
Force[N]
Displacement[mm]
4EL-DC-BD 4EL-FC-BD
Bending dominated Shear dominated
FC DC FC DC
2 EL X X X X
4 EL X X X X
8 EL X (HALF) X(HALF) X X
4.1.7 Load Application
To apply the load to the beam multiple techniques may be used, the most frequently used in
literature were evaluated in this chapter. The evaluated loadtypes are illustrated in Figure 25, a more elaborate figure is presented in APPENDIX VIII.II. The results of the total test population are
presented in APPENDIX IX.VI. These loadtypes were selected to research the influence of 3 parameters:
-Force controlled or Displacement controlled loading -A load applied to the web or to the flange
-A load applied directly to the section or through a reference point that is tied to that section.
FC-DI-DB FC-RN-DB DC-RN-DB FC-RN-DBLOW
FC-DI-FL FC-RN-FL DC-RN-FL
The influence of the loading on the stress distribution was measured. The stress distribution as a result from loading was similar for all models that were loaded through a reference node. Even loads applied directly to the flange or through a reference node resulted in similar stress distributions as illustrated in Figure 26. Even though the stress distribution of the DBlow model in Figure 26 is less smooth, all models resulted in the same failure loads for shear- as well bending dominated analyses.
Figure 26- shear stress distribution over the width of the loaded flange
The models that were loaded at their flanges showed premature local failure for high shear stresses, as is illustrated in Figure 27. Therefore, the flange imposed loadtypes were excluded from further testing.
Figure 27- Premature local failure of flanges due to imposed loads
The concentrated loads combined with the bending moments were scaled proportionally by the same load proportionality factor (LPF). To obtain results at specific V/Vpl ratios, an estimation had to be made for the required loading. During displacement controlled loading it proved hard to predict the rotations and lateral displacements to apply to obtain results for specific M-V combinations. Therefore, force controlled loading was preferred. As was concluded in chapter 4.1.6, the FC models tended to have convergence issues, however, no convergence issues were encountered with FC-DI DB.
-4 -3 -2 -1 0 1 2 3 4 5
0 50 100 150
τ12[N/mm2]
Width of loaded flange [mm]
FL-DI-FC-BD FL-RN-DC-BD DBLOW-RN-DC-BD
In addition to the loadtype, it was examined if the sequence of loading resulted in order effects. Two methods may be used to obtain the bending moment capacity at a certain shear force:
-Increase the concentrated forces and bending moments simultaneously in a single Riks step
-Increase the concentrated forces up to the desired V/Vpl ratio in a static general step, and subsequently increase the bending moments in a Riks step to determine the residual capacity.
Applying the loads subsequently ensured the bending moment capacity was obtained at the exact intended shear force, however, the results could differ from loading simultaneously as nonlinear material behavior, and instability might cause order effects. Figure 28 indicates no order effects were apparent as both methods obtained the same maxima. However, for high shear forces failure can occur prematurely in the static general step. Subsequently the second step was not initiated and the beam was loaded proportionally nonetheless as is illustrated in Figure 28 . Therefore, all simulations were loaded simultaneously to be consistent.
Figure 28- Comparison of 4-point bending test loaded subsequently(dashed) and simultaneously(continuous) 0
0,2 0,4 0,6 0,8 1
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9
M/Mpl[-]
V/Vpl[-]
4.1.8 Imperfections
According to EN-1993-1-5 ANNEX C both geometric and structural imperfections should be included in the FE-model. Imperfections may either be obtained from a refined buckling analysis, or from equivalent values. Any type of imperfection should be taken such that the lowest resistance is obtained. In this section various types of imperfections and methods were evaluated
Geometric imperfection
A geometrical imperfection consists of a combination of an imperfection shape and amplitude. The shape can be equal to one or multiple eigenmodes, or an equivalent shape according to Figure C.1. EN-1993-1-5 ANNEX C as illustrated in Figure 29. Both methods are used to provide node coordinates of the deflected shape, that are subsequently superimposed on the node coordinates of the GMNIA.
(a) (b)
Figure 29- Imperfections for flange(a) and web(b) according to Table C.2. EN-1993-1-5 ANNEX C
To obtain the eigenmodes of the beam a Linear Buckling Analysis(LBA) with a subspace eigensolver was used. Even though this method is widely adopted, a problem with the LBA was encountered for models loaded in M-V interaction. The graph in Figure 30 indicates an increased bending moment capacity for V/Vpl=0,4 to V/Vpl=0.6. Further inspection into these analyses indicated a discrepancy in the shape between the first 14 eigenmodes and the failure mode in the GMNIA, as illustrated in Figure 32 and Figure 31.
Figure 30- M-V interaction graph with imperfections obtained from a LBA λf =10 λw =83 S235 0
0,2 0,4 0,6 0,8 1 1,2
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8
M/Mpl[-]
V/Vpl[-]
Figure 31- First eigenmode of λf =10 λw =83 S235 V/Vpl=0.4 Figure 32- Failure mode of λf =10 λw =83 S235 V/Vpl=0.4
The imperfection obtained from the LBA was not located at the point of failure in the GMNIA; the mismatch resulted in the increased bending moment in Figure 30.
The cause of the mismatch was found when the displacements over the increments of the GMNIA was considered in Figure 33. When stresses are elastic the displacements correspond to the eigenmode of the LBA, subsequently, when the compression flange starts to yield the flange fails premature to buckling of the web. The assumption the member fails in plasticity or elastic buckling seems to be incomplete; elasto-plastic buckling could occur and seems to be sensitive to the applied imperfections.
Obviously, non-linear material behavior is not incorporated in the LBA, and therefore it might provide imperfections in other sections where failure would occur with non-linear material behavior. To check whether the eigenmode is appropriate as input for the imperfection the eigenvalue should be below the elastic moment resistance, however, if not, the researcher cannot switch between methods in one test population.
Increment 9 Mises:
167 N/mm2
Increment 27 Mises:
240 N/mm2 Increment 70 Mises:
242 N/mm2
Figure 33- Evolution of displacements U1 over the increments of the GMNIA
Imperfections should be kept constant to prove causality between slenderness, the influence of shear stress, and moment capacity. This poses another challenge when using eigenmodes as input for the GMNIA analyses. The Abaqus manual recommends applying multiple eigenmodes when the corresponding eigenvalues are close together, however, if the eigenmodes overlap on a single node the imperfections are added up. Therefore, to maintain a constant amplitude the multiple eigenmodes should be selected very carefully. Nonetheless, a constant amplitude for the imperfection for every model is nearly impossible. Even though the LBA is used widely, the difference between linear and non-linear behavior makes it an unsuitable method for many GMNIA analyses.
To obtain the node coordinates of the equivalent shapes a static general step is used. Imposed displacements equal to the amplitude illustrated in Figure 29 were applied to points on the flange and web of the beam.
2 Types of imperfections are applied to the sections in M-V interaction: one shape corresponding to local buckling. one shape corresponding to shear buckling. The stress distribution in a 4-point bending test is identical in both sections, therefore, for any M-V interaction the beam can fail in the lowest failure mode; shear or axial buckling. The amplitude of the imperfection in the outer sections is set to 95% of the middle imperfection, to induce local buckling in the middle and prevent erratic behavior for pure bending.
The buckling shape for shear- and local buckling were recreated from observed failure modes from multiple GMNIA analysis. The double sine is used as buckling shape as it resulted in the lowest bending moment capacity in the GMNIA, as is showed in Table 6.
Table 6- Bending moment capacity for single and double sine imperfection obtained from λf=10 λw=83-S235
Single sine 0.96 M/Mpl
Double sine 0.93 M/MplF
Double sine 0.93 M/MplF