Application of the generalized nonlinear constitutive law in 2D shear flexible beam structures

The paper presents a modified finite element method for nonlinear analysis of 2D beam structures. To take into account the influence of the shear flexibility, a Timoshenko beam element was adopted. The algorithm proposed enables using complex material laws without the need of implementing advanced constitutive models in finite element routines. The method is easy to implement in commonly available CAE software for linear analysis of beam structures. It allows to extend the functionality of these programs with material nonlinearities. By using the structure deformations, computed from the nodal displacements, and the presented here generalized nonlinear constitutive law, it is possible to iteratively reduce the bending, tensile and shear stiffnesses of the structures. By applying a beam model with a multi layered cross-section and generalized stresses and strains to obtain a representative constitutive law, it is easy to model not only the complex multi-material cross-sections, but also the advanced nonlinear constitutive laws (e.g. material softening in tension). The proposed method was implemented in the MATLAB environment, its performance was shown on the several numerical examples. The cross-sections such us a steel I-beam and a steel I-beam with a concrete encasement for different slenderness ratios were considered here. To verify the accuracy of the computations, all results are compared with the ones received from a commercial CAE software. The comparison reveals a good correlation between the reference model and the method proposed.


Introduction
For decades a finite element (FE) analysis has been a popular method for modelling advanced engineering problems. The FE models in comparison to the analytical ones have wider applicability and universality, therefore it is often implemented in various modern engineering tools. Most civil structures can be analysed using simple 2D beam or frame models without sacrificing the accuracy of the results obtained with these simplified models. Therefore, majority of commercial software for strength analyses of engineering structures use the beam, truss or frame finite elements formulation.
Nowadays, the developers of FE software more and more often provide their users new functions extending its capabilities. Following this trend, the software often allows to include users' material or element subroutines tailored for particular needs. However, despite its popularity, they usually allow only linear elastic analysis. In order to extend its functionality with nonlinear (material or geometric) analyses through user subroutines the computational and / or material mechanics expertise is required. A practical solution may be the use of classic beams and frames FE linear solver extended with a simple implemented generalized nonlinear constitutive law (GNCL) algorithm. This idea was first time mentioned by Mahin and Bertero in 1977 [13]. The idea of dividing the cross-section into layers was then used to analyse the strength of reinforced concrete columns. A similar approach was used by El-Tawil and Mirza [3,14], where e.g. uni-and biaxial bending strengths of composite short columns were analysed. In 1978, the layer model was used by Rotter and Ansourian to analyse the behaviour of bending composite beams [18]. The theoretical values were compared with the results obtained from the experiments, obtaining a good correlation. In 1982, Łodygowski used generalized nonlinear constitutive law (GNCL) method for geometrically and physically nonlinear analysis of beams and plane frames [11]. Later, Łodygowski and Szumigała used the division into layers in a two-stage analysis of bending composite beams [12]. In the first stage of the method, the cross-section is discretized and the constitutive law is formulated in the form of bending moment-curvature relationship. In the second stage, the constitutive law is adopted in the finite element nonlinear code.
The two-stage approach was also used by Szumigała [21] to analyse composite steel-concrete frame structures. In this case, the constitutive law was formulated in the form of bending stiffness-curvature relationship. In 2019, Grzeszykowski and Szmigiera used the GNCL method to compute the nonlinear longitudinal shear distribution in composite steel-concrete beams [8].
The proposed method can have many applications in engineering computations regarding multi material cross-sections. It can be used for strength analysis of engineering structures. An example of such application can be found in the paper of Farhan [4], where behaviour of concrete-filled steel tube composite beams was experimentally tested. Also, the proposed method can be adopted as a method of homogenization in cases, in which a cross-section composed of several materials is used.
One of such methods was proposed by Siwiński [19], where the cross-sections of reinforced concrete elements were homogenized. Notice that when using composite structures, it is important to ensure homogeneity in transferring the loads by material integrity. This problem was addressed by Jayanthi [10] where the performance of different types of shear connectors in steel-concrete composite construction was analysed. The GNCL method can be adopted in such cases, but also in all kinds of beam or truss structure analyses. One of the examples is the work by Barszcz [2], in which a multistorey steel structure is analysed. Having in mind that for some complex structures or materials, as examples discussed above, it is important to take into account the shear effect in the finite element model, the paper presents the following method that takes into account the shear forces.

Normal and shear strains
The proposed method is embedded in the classical framework of FE analysis, see Classically, the global stiffness matrix is assembled by considering the stiffness matrices of all elements. In the method proposed, the element stiffness is iteratively reduced during the analysis due to deformations, namely, 0 − normal strains, − shear strains and − curvature, which are computed from the nodal displacements , see Fig. 1a. Normal strain, 0 is taken as the ratio of an elongation Δ to a beam length , which takes the following form: where 1 , 2 are the nodal displacements along the beam axis.
Shear strains, , according to Timoshenko theory are computed as the difference between the nodal rotation and the first derivative of vertical deflection, , of the beam: In the paper, for comparison only, Bernoulli beams are also considered. Bernoulli's hypothesis assumes that the cross-section is perpendicular to the axis of the deformed beam, consequently, this assumes that the shear strains, , are equal to zero. The curvature, , is derived from the beam deflection function ( ), which is represented by a third degree polynomial: The constants of the polynomial are determined from the boundary conditions. Thus, let assume that: From above, we obtain the following: Knowing the beam deflection function, ( ), enables determining the curvature. In the case of small displacements, the following simplification can be taken: Therefore, the curvature reduces to: In the approach proposed here, the representative curvature is used, namely, a weighted mean curvature, ̅ , computed in three Gauss points: where 1 , 2 , 3 are weights, which takes the values: 1 = 5/18, 2 = 4/9 and 3 = 5/18. The curvatures 1 , 2 , 3 are computed in the following Gaussian coordinates: 1 = 1/9 , 2 = 1/2 , and 3 = 8/9 .

Stiffness reduction
In Fig. 1, an overall algorithm of the method proposed here in the form of a block diagram is presented. The method proposed, as an extension of the approach proposed by Szumigała [21], is introduced in the following section. The extension is including the influence of shear strains on a stiffness reduction. Similar to the original version of the method, here, the element cross-section is also divided into thin horizontal layers, see where: , and are normal strains and , and are shear strains. Later, the reduced stress, , in each layer is determined from the reduced strain based on vs. plot for a particular material. If multi-material cross-section is considered, the computations are performed for each material.
In the next step, see Fig. 1b, the Young's modulus, , is computed from the values of stresses and strains by the following: Further, the shear modulus for isotropic materials, , can be computed from: where is Poisson's ratio. After determining the Young's and shear modulus, it is possible to compute the tensile and shear stiffnesses, and , respectively: where and are number of layers in the cross section and number of materials, respectively; and are the total number of layers and number of materials, respectively. is a cross-section area of -th layer and -th material, and ̅ is a shear correction factor for -th material, which is computed from: where is a moment of inertia about the horizontal axis of the cross-section, is a static moment about the horizonal axis of the severed part, is a cross-section area and is a width of a layer. A position of neutral axis in -th iteration of computing is obtained from: Knowing the position of the neutral axis enables computing a moment of inertia of -th layer, ; then a bending stiffness, , may be determined according to the expression: In the original version of the method [21], a cross-section analysis is done prior to main computations.

Materials
The proposed method allows to use a nonlinear constitutive law of any material.
In the examples analysed here, the nonlinear law describing the behaviour of concrete and steel were used.
In Table 1, the engineering parameters of the materials used in the study are presented, where is a medium value of a compressive strength of concrete, is a characteristic value of a tension strength of concrete, is the fracture energy and is a yield strength of steel. In Fig. 2, the nonlinear stress-strain relations for concrete and steel are presented. The reference FE models used in the results section were computed in a commercial software of ABAQUS FEA [1].
In these examples, the steel was described by an elastic perfectly-plastic model: where is the Young's modulus of steel.
For concrete in compression, a nonlinear model presented in Eurocode 2 was used: where = / 1 and * is computed as (2.16) where is the Young's modulus of concrete and 1 is given in Eurocode 2 for a different grades of concrete.
In order to simulate behaviour of concrete in compression and in tension, a concrete damage plasticity (CDP) model with additional definition of a fracture energy was used: where ̅ is the effective stress and * is a scalar degradation parameter. More formal details and practical numerical outlines regarding CDP model may be found in the paper of Jankowiak and Łodygowski [9]. An identification procedure of main concrete parameters for this constitutive model can be found in Gajewski and Garbowski [5], alternative models can be characterized using approach described in the paper of Gajewski and Garbowski [6], Garbowski et al. [7] or Zirpoli et al. [22]. The material parameters used in the study for the concrete in CDP model are specified in Table 2, while properties of the steel are already presented in Table 1.

The analysis of bending, tensile and shear stiffness
Apart the possibility of computing the internal forces of particular beam structure by using GNCL, it is also possible to determine stiffness reduction plots for defined cross-sections and adopted material models prior to FE computations. The stiffness reduction for three examples of different crosssections is presented here, namely, steel IPE240, encased IPE240 and encased IPE450 with a concrete slab.
In Fig. 3, the selected diagrams of bending, tensile and shear stiffnesses, , and , for the IPE240 steel cross-section depending on the internal forces are shown. In Fig. 3a, it may be observed that the stiffness is much more sensitive to the change of the bending moment, , than to the change of the normal force, . Fig. 3c shows that stiffness has comparable sensitivity to the change of the shear force, , and the bending moment. In all cases, the double symmetry is observed.

Bernoulli vs. Timoshenko beam
In order to verify the effect of taking the shear into account on the beam displacement the relations between results from Timoshenko and Bernoulli theories were compared. In Fig. 6  Timoshenko beam displacements (new approach) to Bernoulli beam displacements (classic approach [21]) is shown. In results, as expected, the smaller the slenderness ratio of the beam, the greater the influence of the shear force on the displacements. For slenderness ratio of 20, the effect of the shear force is less than 1% and increases hyperbolically to almost 10% for a slenderness ratio of 5.

Proposed method: Timoshenko theory vs. Bernoulli theory
To demonstrate the influence of the shear force on the behaviour of the structure, 1.2 m, 1.8 m, and 2.4 m long beams were analysed here, including (with Timoshenko theory) and excluding (with Bernoulli theory) shear effect in the GNCL model. All computations in this subsection were computed using GNCL method and by using Model 1 (see Section 2.5).
In Fig. 7, the diagrams of the force in the mid-span of the IPE240 beam with respect to the displacement applied are shown. The greater the length of the beam (more slender structure), the smaller the difference between the forces obtained for the Bernoulli and Timoshenko theories; this is due to lower influence of the shear force on the displacements, see Fig. 6. In all cases, as expected, the forces for the Bernoulli beam were greater than those obtained for the Timoshenko beam. In Fig. 7a, force-displacement curves for the beam with the shortest length are shown. In this case, as expected, the differences between the curves obtained by the Bernoulli and Timoshenko theory are the greatest. For the 1.8 m and 2.4 m long beams, shown in Fig. 7b and Fig. 7c, the differences are smaller because the influence of the shear force is smaller for the beams with greater slenderness ratio. The percentage differences between ultimate load decrease with increasing beam length and are equal 8.7 %, 6.4 % and 2.2%, respectively.
In Fig. 8, the force-displacement diagrams in the mid-span of the encased IPE240 beam are shown.
As in the previous example, the forces obtained by Bernoulli theory are greater than for Timoshenko theory. Fig. 8 shows that the shorter the beam, the greater the differences between values obtained by the two theories. The influence of the shear force on the load capacity is smaller and equals 3.  In Fig. 10, the force-displacement diagrams in the mid-span of the encased IPE240 beam are shown.

GNCL method compared to elasto-plastic approach
Notice that, here, the reference model, a 3D solids FE analysis was used (Model 2, see Section 2.5).
After plasticity was reached, the static equilibrium paths stabilized at a certain level. As in the example of the IPE240 beam, it may be observed a greater difference for the length of 1.2 m than for the lengths of 1.8 and 2.4 m. The differences in load capacity are 8.1 %, 4.2 % and 3.9 %, respectively.
The results obtained by the method proposed and FE commercial software are in good agreement (see Fig. 9 and Fig. 10). Also, the classic and proposed GNCL approach, however, for Bernoulli theory only, were compared (not included here); the mean error between forces obtained for encased Moreover, structural modelling is easier because it does not require modelling experience, nor expertise knowledge of advance techniques of the finite element method, which would allow to simulate the structure with complex materials or geometry. In particular, the proposed method can be used in beam and truss structures, such approach eliminates the need to model 3D elements taking into account the contact in cross-sections composed of several materials with nonlinear constitutive laws. The method seems to have many promising applications in civil engineering and mechanics, for instance in modelling structures with shear connectors, multi-layers or sandwich panels.

Conclusions
In the paper, the method of generalized nonlinear constitutive law (GNCL) together with FE formulation including shear theory in beams was presented. It enables computing the shear strains and take them into account in the element stiffness reduction. The GNCL method derived here enables a nonlinear description of materials used, i.e. concrete and steel. The stiffness of various cross sections depending on internal forces was shown. In the presented examples, the beam structures with a different length to height ratios were analysed. Performed computations for various beam slenderness ratios with Timoshenko beam theory showed expected influence of the shear force on the structure behaviour (deflections).
The method proposed allows an easy consideration of material nonlinearities in the beam/frame models. By applying a GNCL model, computations can be performed for complex cross-section composed of several materials with different physical properties. This may be obtained not only for slender structures, but also in cases of short beams, in which the shear effects are crucial. The GNCL provides the simple homogenization of the complex cross-section which then enables the iterative stiffness reduction based on the internal forces or deformations. The static equilibrium path in a nonlinear form can be easily computed with this approach. Main advantage of the method is that there is no need to build a full 3D model with elasto-plastic constitutive models and iterative solvers. This makes the modelling of the structure easy; one does not require a knowledge and experience in FE method. Also, compared to commercial software, the method gives good results and the computational time is several times shorter, due to small number of degrees of freedom, while comparing with 3D models.