FEM modelling of the static behaviour of reinforced concrete beams considering the nonlinear behaviour of the concrete


 This paper presents a finite element (FE) method of modelling reinforced concrete beams. The proposed model takes into account the phenomena characteristic of reinforced concrete structures, such as the interaction between two materials (concrete and steel), the cracking caused by mechanical loads and the variation of the Young's modulus under increasing load. A relevant numerical FE analysis was carried out in the ABAQUS system using the concrete damaged plasticity (CDP) model. The character of Young's modulus variation due to increasing stress intensity level was taken from the author's own research. The results of the FE calculations were compared with the results yielded by the author's numerical bar model.


Introduction
Nowadays structural engineers are often adopting innovative design solutions not covered by the schemes and algorithms proposed by universal design codes. As part of the design process, object models are created and subjected to analyses. Formerly mainly experimental (reduced-or full-scale) models and simple mathematical models were used chiefly due the fact that numerical techniques were then little advanced and not generally available. The rapid development of electronics and informatics has impacted other scientific technical fields, including civil engineering. Owing to this, one can now create and analyse complicated numerical models (especially FEM models) whose cost is relatively low (in comparison with experimental models and the associated measuring apparatus).
In numerical modelling, the adopted assumptionsparticularly the strength and deformation parameters of the materials of which the structure is to be made -are of utmost importance. This paper presents a procedure for the FEM modelling of reinforced concrete beams, using concrete parameters obtained from in-house tests (mainly the Young's modulus). It can be said that this is one of the most effective ways to reduce experimental costs. In the considered case, only material tests were carried out, which is not sufficient, because the FEM model should be confronted with tests of model of the structure in half or full scale. However, FEM modelling can reduce the need for multiple structural models. The created numerical model took into account phenomena characteristic of reinforced concrete structures, such as the interaction between two materials (concrete and rebars), cracking and the variation of the Young's modulus of the concrete depending on the stress intensity level. The variation was recorded in in-house tests.
In other words, the main goal of this study was to implement the experimental results in the created numerical FEM model. Moreover, a reinforced concrete beam bar model, which takes into account the Young's modulus changing depending on the stress intensity level, was developed in order to verify the results. stiffness and the stiffness of elements in a spatial stress state.
Young's modulus applies to materials for which the stress-strain relationship is linear. Although the behaviour of concrete is nonlinear, according to Neville (1977) it can be assigned the following Young's moduli in compression: -tangent initial Young's modulus -defined by a straight line tangent to the graph in the origin of the coordinate system ( Fig. 1 -straight line 1); -tangent Young's modulus -defined by a straight line tangent to the graph in a given graph point ( The slope factors of the straight lines (1,2,3,4) are the Young's moduli corresponding to the respective straight lines. The Young's modulus of concrete is a very complex problem conditional on many factors. According to Neville (1977), these include: -the aggregate (the Young's modulus of the aggregate, the shape of its larger particles and their surface area), as described by Jurowski and Grzeszczyk (2015); -the cement paste/aggregate ratio (the Young's modulus of the aggregate is usually higher than that of the cement paste); -the age of concrete (the older the concrete, the higher its Young's modulus); -the moisture content of the specimen (the Young's modulus of a moist specimen is higher than that of a dry specimen); -the bulk density of the concrete (the higher the bulk density of the concrete, the higher its Young's modulus); -the temperature of the initial curing of the concrete (the Young's modulus is higher when the concrete is subjected to a lower temperature during its initial curing).
Moreover, according to Neville (2011), the effect of the stress intensity level in concrete on the magnitude of tangent and secant Young's moduli is stronger for concretes with lower compressive strengths. Such concretes behave more ductilely. Whereas concretes with higher compressive strengths are more brittle.
In EN 12390-13 (2013), one can find two methods of testing the secant elasticity modulus of concrete in compression, i.e. method A and method B. In both the methods the specimen is preliminarily loaded in a few cycles up to a specified stress intensity level (0.1÷0.15 f cm in method A and 1/3 f cm in method A and B, where f cm is the mean compressive strength of concrete in a uniaxial stress state, determined for collocated samples) and subsequently compressed until failure. The specimen setting in the strength-testing machine is assessed on the basis of strain readings from the relevant test stages, and appropriate adjustments are made to make sure that the eccentricity in compression will be as small as possible (the specimen should be compressed axially). Using the stress and strain obtained from the particular load cycles, one calculates the secant initial Young's modulus (method A) or the secant stabilized Young's modulus (method A and B).
According to ASTM C469/C469M-14 (2014), the specimen should be subjected to at least three load cycles. The purpose of the first load cycle is to check and if need be, to correct any deviations from the proper specimen setting in the strength-testing machine. In the second and subsequent cycles, the stress at which the strain amounts to 0.05‰ and the stress corresponding to 40% of the mean compressive strength of the concrete in collocated samples are read. On the basis of the read stress and strain, one calculates the Young's modulus of the concrete as a mean for the particular cycles.
In the case of concrete, one can also distinguish the so-called dynamic Young's modulus which is determined not through a static compression test, but by investigating longitudinal, transverse (bending) or torsional vibrations, as described by Neville (1977). Musiał and Grosel (2016) tried to investigate the dynamic Young's modulus of concrete by exciting the vibration of concrete beams with different reinforcement ratios, made using different concrete mix designs. According to Neville (1977), this Young's modulus is considered to be approximately equal to the tangent initial Young's modulus determined in the static compression test. This is so since only small stress arises in a vibrating specimen. However, the relation between the static Young's modulus and the dynamic Young's modulus is not a simple one. Jurowski and Grzeszczyk (2018) presented a tentative relation between the two quantities.
Knowledge of the Young's modulus is very important in order to design civil engineering structures. The calculations can be carried out according to the nonlinear analysis, where the stress-strain relationship for the concrete is used. However, where the linear-elastic analysis is usually sufficient, the Young's modulus of the concrete is used. For example, in the case of beam elements, knowing the Young's modulus, the element's dimensions (the dimensions of its cross section, and its length), its load (the axial forces and the bending moments) and the boundary and continuity conditions, one can calculate the element's contraction or elongation under, respectively, compression and tension and its deflection under bending. Moreover, axial stiffness and bending stiffness have a bearing on the distribution of internal forces in statically indeterminate systems, while bending stiffness influences second-order effects which often should be taken into account in the design of structures. Considering the above, the author decided to analyse the static behaviour of reinforced concrete beams, taking into account the Young's modulus of concrete changing depending on the latter's stress intensity level, i.e. taking into account the nonlinear behaviour of the concrete. This paper presents a way of modelling a reinforced concrete beam in ABAQUS using the concrete damaged plasticity (CDP) model which takes into account the nonlinear behaviour of the concrete, the cracking of the reinforced element and the issue described in detail by Pędziwiatr (2008Pędziwiatr ( , 2009) -the interaction between concrete and rebars. Similar problems have already been

Materials
Tests were carried out on 12 cylindrical concrete specimens 300 mm high and 150 mm in diameter. The concrete mix design is presented in Table 1.
Concrete mix consistency was controlled by gradually adding superplasticizer Sikament FM-6. Ultimately, the superplasticizer mass was assumed to be equal to 1.2% of the cement mass.

Experimental procedure
Investigations of the Young's modulus of concrete were carried out using two methods (I and II) based on Code method B (in accordance with EN 12390-13, 2013) for determining the stabilized Young's modulus of concrete. The method was modified to take into account the influence of different concrete stress intensity levels on the Young's modulus.
Three electrical resistance strain gauges RL=350/50 were stuck at every 120° on each cylindrical specimen at half of its height along three lines on the cylinder side wall. The shape of the specimens and the arrangement of the electrical strain gauges are shown in Fig. 2. In the case of method I the load was applied according to the loading diagram shown in Fig. 3.
The symbols used in the description below should be interpreted in accordance with Fig. 3.
The numbers in brackets (1)-(6) represent the particular stops, i.e. pauses no longer than 20 s during which the specified stress level in the specimen should be maintained.
Method I differs from the standard code method in this that mean stress and strain in the specimen in the final load cycle should be read also halfway between lower stress level σ p and upper stress level σ a and above stress level σ a at every 1/6 f cm .
The secant Young's modulus of concrete for a given stress interval is calculated as the ratio of the difference between the stress indicated by two neighbouring readings and the difference between the mean strain corresponding to the readings.
In the case of method II, the load was applied in accordance with the diagram shown in Fig. 4.
The symbols in the description below should be interpreted consistently with Fig. 4. The numbers in brackets (1)-(12) represent the particular pauses.
In the case of method II, in the first two cycles the specimen is loaded up to the stress level of 2/3 f cm at a stop at 1/3 f cm . In the final load cycle, the mean strain and stress in the specimen are read at the end of each stop. The readings should be taken for stress halfway between stress levels σ p and σ a and between stress levels σ a and σ a ', whereas after stress σ a ' is reached and a stop is made the specimen is loaded at stops at every 1/6 f cm .
In the case of method II, the secant elasticity modulus of concrete for a given stress range is calculated similarly as in method I.
The subsequent photographs show the specimen in the strength-testing machine prior to the test (Fig. 5) and after failure (Fig. 6). The visible characteristic cones in the  core of the failed specimen indicate that the specimen was compressed axially.

Test results
The test results for an exemplary specimen tested according to method I are presented in Figs. 7-9. The specimen loading diagram is shown in Fig. 7. The stressstrain curve drawn on the basis of the test results is shown in Fig. 8. The dependence between the Young's modulus of the concrete and the compressive stress and the stress intensity level in the specimen is illustrated in Fig. 9.
In accordance with EN 1992-1-1 (2004), the secant elastic modulus of concrete for the specimen represented by Figs. 7-9 amounts to E cm =30.28 GPa.   The test results for the exemplary specimen tested according to method II are presented in Figs. 10-12. The specimen loading diagram is shown in Fig. 10. The stressstrain curve is shown in Fig. 11. A diagram of the Young's modulus of concrete versus the compressive stress and stress intensity level in the specimen is shown in Fig. 12.
It appears that method II yields an almost constant value of the Young's modulus for the stress intensity level of about 60%. Whereas the value of the Young's modulus investigated using method I begins to decrease at lower stress intensity levels. According to Neville (1977), the multiple loading of concrete eliminates creep (which occurs already as the concrete is being loaded). Subsequently, the value of the Young's modulus stabilizes and the curvature of the stress-strain curve decreases.  This can be observed for testing according to method II, where in the first two cycles the concrete is loaded not to the stress level of 1/3 f cm as in method I, but to 2/3 f cm . As a result, the material adapts to higher stress intensity levels, as shown in Fig. 12.
Thus for the same concrete not subjected to loading yet versus already loaded one gets different elastic modulus curves. Therefore in a precise analysis it is important whether the considered element is a debuting element or whether it has already been subjected to loading.

Numerical analysis 4.1 FEM model
An analysis of the static behaviour of the beam in ABAQUS was carried out for the concrete tested using method I (the specimen represented by Figs. 7-9). Reinforcement ratio ρ=1% (commonly used and economically viable) was assumed for the beam.
The geometry of the considered element is shown in Fig. 13.
The following nominal loading steps with the load uniformly distributed along the whole length of the beam were considered: 40, 45, 50, 55, 60 and 65 kN/m.
The beam was modelled in ABAQUS, taking into account the CDP model. The stress-strain relationship obtained from the tests was entered into program. Moreover the actual arrangement of reinforcement was taken into account in the model.
The beam alone was defined as a three-dimensional (3D) solid (element type: solid), while the reinforcement was built of one-dimensional elements in 3D space  (element type: wire). Moreover, steel support pads (element type: solid) were designed in order for the FE analysis to map well the actual support.
The rebars and the steel support pads were assigned linearly elastic (elastic) materials. The rebars were assigned reinforcement steel with E = 200 GPa and Poisson's ratio ν = 0.3, while the steel support pads were assigned structural steel with E = 210 GPa and ν = 0.3.
A definition of a material such as concrete, to which a nonlinear stress-strain curve is to be assigned and which is to take element cracking into account, is more complex and requires a special material model, i.e. the CDP model. Moreover, one should endow the concrete with the elastic material characteristic since the stress-strain relationship in CDP is a polyline and the Young's modulus assigned to the elastic characteristic corresponds to the slope of the first (counting from the origin of the coordinate system) segment of the stress-strain polyline. In addition, the elastic characteristic includes a Poisson's ratio which for the concrete was assumed to amount to ν = 0.2.
The following concrete parameters were used in CDP: where: ε c -strain due to compression for a given point on the stress-strain curve, ε 0c el -a part of the strain due to compression, normalized to the initial Young's modulus, expressed by the eq.
where: σ c -compression stress for the given point on the stress-strain curve, E 0 -the initial Young's modulus. The notations in eqs. (1) and (2) should be interpreted in accordance with Fig. 14.
When characterizing concrete in tension, taking into account tension stiffening, one should specify stress and the corresponding cracking strain values. The stressstrain curve for concrete in tension is shown in Fig. 15.
The tensile strength of concrete was assumed on the basis of the following eqs.   Cracking strain is calculated from eq. (5): 0t el 0t where: ε t -strain due to tension for a given point on the stress-strain curve, ε 0t el -a part of the strain due to tension, normalized to the initial Young's modulus, expressed by the eq. (6): where: σ t -tensile stress for a given point on the stressstrain curve, E 0 -the initial Young's modulus.  Moreover, data relating to the failure parameter and the corresponding cracking strain were entered in order to obtain results for the failure of concrete in tension. The failure parameter can be estimated using the eq. (7): where: σ t0 -the tensile strength of concrete in a uniaxial stress state. The notations in eqs. (5)-(7) should be interpreted in accordance with Fig. 15.
Taking into account the above, the stress-strain characteristic (based on the author's own research) shown in Fig. 16 was entered into the program.
Material (bulk) density was assigned to both the steel and the concrete since in the considered case a dynamic explicit analysis which requires this parameter is carried out.
The object model consists of the beam's concrete, the rebars and the steel support pads. Beam-rebars interactions and beam concrete-steel pads interactions were introduced in the next step. The rebars were connected with the beam's concrete through an embedded region constraint while steel pads-concrete interactions were effected through a tie constraint. The model is shown in Fig. 17.
Then boundary conditions and loads are entered. The supports were introduced as linear pinned supports with the possibility of shifting one beam end along the beam. The supports were located at half the width of the steel pads. The load was defined as a pressure-type load applied to the beam's top surface.
Subsequently, all the model's parts were divided into FEs. C3D8R FEs (8-node reduced integration-type cuboid FEs) were assigned to the beam and the steel pads. The element size amounted to approximately 50 mm for the beam and to 25 mm for the steel pads. T3D2 FEs (two-node truss elements in 3D space) with the element size of 50 mm were assumed for the rebars.
After computations, the results for vertical displacements (an exemplary map is shown in Fig. 18) and the images of cracking (see Fig. 19) were analysed.

Numerical bar model
The deflections of the beam analysed in ABAQUS were also calculated using a numerical bar model developed by the author on the basis of the Euler-Bernoulli beam theory. The procedure was carried out for the same concrete which had been assigned to the beam in the FE analysis. Figure 20 shows the beam structure and loading diagram.
The proposed method requires knowledge of beam deflections according to EN 1992-1-1 (2004). For example, the considered beam's deflection calculated in accordance with EN 1992-1-1 (2004) for the load of 40 kN/m amounts to 10.20 mm.
In order to calculate the beam's deflection, taking into account the changes in the Young's modulus, an intermediate cross section (between cross section behaviour in phases I and II) was determined. Such a cross section with a certain concrete tension zone is determined to take into account tension stiffening. When the obtained moment of inertia is substituted into the formula for the deflection of a homogenous prismatic bar, one gets the deflection value calculated according to EN 1992-1-1 (2004). Then the beam was divided into n=40 equal fragments. The bending moment along the length of each of them is constant and equal to the calculated bending moment in the middle of each of the fragments. Consequently, the diagram of bending moments for the whole beam has a stepped shape as shown in Fig. 21.
Each of the fragments has an identical moment of inertia (for the cross section mentioned above), but the Young's modulus for each of them was calculated independently on the basis of the compressive stress in the cross section. The Young's modulus for a given fragment is an integral mean of the Young's modulus of the concrete in the compressed part of the cross section of the considered fragment. The graphical interpretation of the determination of the Young's modulus for the given fragment of the beam is shown in Fig. 22., where: xthe extent of the compression zone, t -the extent of the  tension zone in concrete, ε cg -normal strain in top concrete fibres, ε cd -normal strain in bottom concrete fibres, ε s1normal strain in bottom reinforcement steel, ε s2 -normal strain in top reinforcement steel, σ c,g -normal stress in top concrete fibres, σ c,d -normal stress in bottom concrete fibres, σ s1 -normal stress in bottom reinforcement steel, σ s2 -normal stress in top reinforcement steel, α e -a ratio of the Young's modulus of the steel to the secant Young's modulus of the concrete, acc. to EN 1992-1-1 (2004), E ithe Young's modulus for the i-th fragment of the beam, E 1 -the Young's modulus of the concrete along segment x 1 , E 2,i -the Young's modulus of the concrete at the end of segment x 2 in the i-th fragment of the beam, x 1 -the extent of the compression zone in which the Young's modulus of concrete is constant, x 2 -the extent of the compression zone in which the Young's modulus of concrete is variable.
The Young's modulus of the concrete for a given fragment of the beam was calculated from the eq. (8): 0t el 0t After moduli E i for each fragment of the beam had been calculated, the deflection was calculated by adding up the displacements of the ends of all the fragments in one half of the beam. The displacements were calculated as for a cantilever loaded on its free end with a point moment (the bending moment is constant along the length of the fragment), taking into account the angle of rotation of the cantilever beginning (the cantilever support). For each fragment this angle is equal to the angle of rotation of the end of the preceding fragment. The starting point for the calculations is the first beam fragment (numbering as in Fig. 21) for which the angle of rotation of its beginning is equal to zero (φ 0,1 =0). It is the consequence that the angle of rotation of the axis of a simply supported beam at half of its span for a load uniformly distributed along the beam's length is equal to zero. The formula for the deflection of the beam with Young's modulus variation taken into account is eq. (9): 0t el 0t 0.
The graphical interpretation of the determination of the beam's deflection is shown in Fig. 23. The deflection of the considered beam with the variable Young's modulus of the concrete taken into account amounted to α'=10.26 mm. Similar calculations were made for the loads of 45, 50, 55, 60 and 65 kN/m. The results are compiled in Table 2. Table 3 shows the deflection values calculated using the bar model (taking into account, respectively, the constant Young's modulus of concrete according to EN 1992-1-1 (2004) and the variable Young's model according to the author's own numerical model) and the ones obtained using the ABAQUS software for numerical FE analyses.

Discussion of results
In addition, the data presented in Table 3 are visualized in Fig. 24.
Moreover, the effect of taking into account Young's modulus variability on the design of reinforced concrete beams with regard to the allowable element deflection where: I e -the effective moment of inertia, I g -the moment of inertia of an uncracked cross section, neglecting the reinforcement, I cr -the moment of inertia of a fully cracked cross section, M cr -the cracking moment, M a -the maximum unfactored bending moment in the element. After appropriate transformations of the formulas for beam deflection according to EN 1992-1-1 (2004) the eq. (12), analogous to eq. (11), was obtained:

Conclusions
The results yielded by ABAQUS differ by 2.11-14.88% from the ones obtained using the proposed numerical method for the analysis of a bar beam model taking into account the variable Young's modulus. The higher the load value, the larger the difference (except for loads between approximately 40 and 43 kN/m). This is due to the fact that the presented bar model analysis is an approximation and so it does not fully depict the actual nonlinear behaviour of the concrete in the cross section of a reinforced concrete beam. Moreover, in ABAQUS the beam was modelled as a 3D solid, whereby the deflections yielded by the FE analysis include deformations generated by all the stress components occurring in the element, whereas the proposed method is based on the Euler-Bernoulli beam theory and so takes into account only deformations generated by stress component σ x .  Moreover, it has been shown that when the variability of the Young's modulus is taken into account, this has no significant effect on the design of reinforced concrete beams as regards the allowable deflection since when this variability is taken into account, the relative increment in deflection amounts to about 1.5% (in the case of the considered element). This applies to both the calculations done according to EN 1992-1-1 (2004) and the ones done according to ACI Code 318-19 (2019) since in the analysed case the two procedural algorithms yield similar results.
The aim of this research was to combine classical experimentation with numerical FE analyses. Tests were carried out on relatively small and inexpensive cylindrical concrete specimens and the concrete's mechanical parameters determined in this way were used to create a 3D reinforced concrete beam model which took into account the nonlinear behaviour and failure of the concrete. As a result, it became possible to carry out simulations of the behaviour of the considered beam under mechanical load. Thanks to this approach, expenditures entailed by carrying out tests on the real element were avoided. Material tests are not enough. It is important to calibrate the FEM model on the basis of tests of model of the structure in half or full scale. The results of the FEM analysis can then be considered as reliable. Therefore, it should be remembered that the presented analysis is an introduction to further considerations. In addition, in the further beams' tests the uniform load should be replaced with three-point or four-point bending.