Selection of appropriate concrete model in numerical calculation

Modelling of composite and reinforced concrete structures requires very precise determination of material parameters and constitutive relations between strain and stress. Erroneous selection of the dependencies and incorrect modelling, in particular, of the performance of concrete in tension may generate results in finite element method programs, that are far from the results obtained in an experiment. Using the example of a concrete damage plasticity model, based on fracture mechanics, this paper describes the physical interpretation and the method of the selection of parameters necessary for the appropriate modelling of concrete in a complex stress state. The correctness of the assumed description of concrete was verified on the basis of results of laboratory tests. A comparative analysis of the experimental and numerical results showed that the application of the concrete damage plasticity model allowed correct determination of the concrete element damage mechanisms for each level of strain.


Introduction
The appropriate modelling of composite structures [1] and reinforced concrete structures [2] is a major challenge, in particular due to the material properties of concrete.Unlike steel, concrete is a material that -in compression -behaves in a non-linear way from the very beginning.In addition, it undergoes degradation in a much faster manner [3] when in tension.Hence, the correct description of the concrete failure is the primary issue when considering the behaviour of concrete in a structure with advanced strains.As a result of the development of existing micro-cracks in concrete, the material softens.On the grounds of mechanics of continuous mediums, this phenomenon is defined as damage during which the tangent stiffness matrix ceases to be defined positively.In the case of concrete tension, due to the development of the damage, the tangent stiffness matrix is negatively defined, which leads to the location of the strains.From the mathematical point of view, a change occurs in a type of partial differential equations, which govern the process, from elliptic to hyperbolic.The problem is ill-conditioned, and consequently, the obtained solution depends on the density of the applied discretisation [4].
The presented problem can be regularised via different methods, among which one can distinguish two main groups.The first group concerns the regularisation at the level of a mathematical formulation, while the other regularises the problem at the level of a numerical formulation.The first group includes, above all, viscoplastic models in which the stress state depends on the rate of the strain [5].They are used mainly in dynamics.This group also includes the Cosserat and micro-polar models which are used for soils and granular media [6].In the case of modelling a failure in concrete, non-local models are also used [7], as well as the socalled higher-order, or gradient, models [8].The nonlocal models introduce an averaging function which converts local variables into non-local ones, according to specific weights.
In turn, the possibility of regularisation at the numeric formulation level consists in the explicit introduction of the width of the plastic strain location zone into a finite element [4].In addition, an automatic refinement of finite element mesh is used, based on local errors caused by large gradients of internal variables, or a displacement field discontinuity within the finite element is introduced [9].This leads to the regularisation of the problem.All the above regularisation methods bring some characteristic internal length scale, which provides information on the width of the strain location zone in the material model.

Concrete damage plasticity model 2.1 Mathematical description
The numerical analysis of the initial boundary value problem with the location of failure requires a comprehensive, constitutive modelling of the material.In this model, the scalar damage parameter is used for modelling of the failure of concrete separately in compression and tension.The scalar description of the material damage [10] was initially used to describe the phenomenon of creep.The primary issue in the numerical analysis of concrete structures is to define the correct mechanism of structure failure and its load carrying capacity.This requires the identification of constitutive parameters of the model based on specific laboratory experiments [11].
The constitutive equation [10] of the material with scalar damage is as follows: where: σ -the Cauchy stress tensor, d -scalar damage parameter (representing stiffness degradation), ε -strain tensor,  0 el -initial tensor of elastic constitutive stiffness,  el -degraded tensor of elastic constitutive stiffness.
It is necessary to specify an effective stress tensor based on the equation ( 2): where: ε plplastic strain tensor, and the evolution of the damage parameter d as follows: This, variable d is assumed to be a function of the effective stress  and the equivalent plastic strain  pl .In the concrete damage plasticity model, the stiffness of the material is initially isotropic, but during the process of stiffness degradation, it is determined by two variables: d c in the compression and d t in the tension.Ultimately, the Cauchy stress tensor σ is proportional to effective stress tensor , and the coefficient of proportionality is (1 -d): The damage of the material in compression and tension is determined independently by two variables,  c pl and  t pl , respectively, which specify the equivalent strains (in compression and tension): Cracking (in tension) and crushing (in compression) in concrete are governed by increasing hardening (softening) variables.These quantities control the evolution of the loading surface F and the material stiffness degradation.
The loading surface determines the strain state as a function of stress and effective plastic strain.In the case of non-viscous damage plasticity model, the state of stress and strain must satisfy the condition: Plastic flow is determined by the plastic potential function G(σ) and a non-associated flow rule in the form:

Strength hypothesis and its parameters
The failure of concrete can be divided into two types.The first type refers to cracking in tension, which is associated with the occurrence of a dominant crack and the decline in load-carrying capacity in the direction perpendicular to it [12], while the other type is related to crushing in compression, during which many small tight cracks develop and evolve as a result of the loss of most of the concrete strength related to spalling.The strength of concrete in a complex stress state is dependent on the complexity of the stress state present in the material.Thus, it cannot be determined only based on simple tests, such as uniaxial compression or tension.The strength of concrete can be determined only by considering the interaction of various stress state components.
Therefore, if one desires to describe the strength of concrete with a triaxial state equation, one must present its plane in a three-dimensional stress space, where on the surface there are stress states corresponding to the failure of the material, while inside there are safe performance states.In addition, inside the space there is the so-called plastic potential surface, after exceedance of which two situations occur [13]: • an increase in strains without changing the stress (ideal plasticity), • material softening -destruction of the material.
One of the most widely used hypotheses regarding the strength of concrete is the Drucker-Prager hypothesis [14], under which the failure depends on the shear strain energy and the failure surface in the stress space takes the shape of a cone (Fig. 1 [15]).An advantage of using such a failure criterion is the surface smoothness and thus no complications with numerical applications, but a disadvantage is the incomplete compliance with the real behaviour of concrete.The CDP model used in Abaqus software is a modification of Drucker-Prager strength hypothesis.In the recent years, subsequent modifications have been made in it [16,17] according to which the failure surface in the deviatoric cross-section does not need to be a circle and is regulated by parameter K c .The physical interpretation of parameter K c is the ratio of the distances from the hydrostatic axis to the compression and tension meridians in the deviatoric cross-section.This ratio is always larger than 0.5.If its value is 1, the deviatoric cross-section of the failure surface becomes a circle, as is the case in the classic Drucker-Prager hypothesis (Fig. 2 [15]).
In the paper [13], in accordance with the results of the experimental studies, it is observed that this value is 0.6 and slowly increases with the decrease in tension in the case of normal average stress of equal to zero.In accordance with [15], the recommended value K c is 2/3.This shape is comparable to the strength criterion formulated by William and Warnke [18], which combined three mutually tangent ellipses.It is a theoretical-experimental criterion built based on the results of tests in a triaxial stress state.A similar change applies to the shape of the meridians of the surface in the stress space, which are curved lines.In the CDP concrete model, the plastic potential surface in the meridian plane takes the form of a hyperbole (Fig. 2 [15]).The adjustment of the shape is done via the eccentricity plastic potential parameter.This is the length of the segment measured along the hydrostatic axis between the hyperbole vertex and the intersection of the asymptotes of this hyperbole (the centre of the hyperbole).The eccentricity parameter can be calculated as the ratio of the tensile strength to the compressive strength [19].For the CDP model, the recommended value is є = 0.1.In a limit case where it is equal to 0, the surface becomes a straight line in the meridian plane -the classic Drucker-Prager hypothesis.[15].

Fig. 3. Strength of concrete under biaxial stress in CDP
Another parameter that describes the behaviour of concrete is the point at which it fails in biaxial compression (Fig. 3 [15]).The value σ b0 /σ c0 (f b0 /f c0 ) determines the ratio of the strength of concrete in a biaxial state to the strength in the uniaxial state.
The most reliable research in this area is the study carried out by Kupfer in 1969.After its approximation with an equation of an ellipse, the uniform biaxial compressive strength f cc = 1.16248 [13] or by default 1.16 according to [15].
The last parameter describing the behaviour of concrete in a complex stress state is the dilation angle, that is the angle of the inclination of the asymptote of the failure surface in relation to the hydrostatic axis, measured in the meridian plane (the angle of internal friction in the concrete).Numerical analyses usually assume the value of ψ = 36° or ψ = 40°.

The stress-strain curve in uniaxial compression and tension
The most accurate method to describe the relationship between stress and strain is to perform uniaxial compression testing for the applied concrete and then transform the variables from the obtained graph.In the CDP model, one must determine the so-called inelastic strain  c in by deducting the elastic portion that corresponds to the undestroyed material from the total strain (recorded as part of a uniaxial compression test): During the transformation of strains, one must determine the moment from which the material is defined as non-linearly elastic (Fig. 4 [15]).[15].

Fig. 4. The stress-strain curves of concrete in uniaxial compression
Laboratory tests -uniaxial tests -show that such behaviour of material occurs almost from the beginning of the process of compression, but it is negligible in the initial stage for most numerical analyses.In the paper [13], the author writes that the end of linear elasticity should grow with the increase in concrete strength and should be more accepted from there, rather than determined on the basis of tests.In addition, the author described it as a percentage of stress in relation to the strength of concrete, according to the relationship: A simpler method for determining this ceiling is to arbitrarily adopt it as 0.4•f cm .In accordance with [20], the modulus of elasticity of concrete is defined as secant within the range of 0-0.4•f cm .It is good at this stage to adopt such a ceiling of the beginning of the inelastic phase so that the initial value of the Young's modulus and the secant value are consistent with each other.Most numerical analyses usually examine the phase of reaching the strength of the material, rather than the initial phase of its performance.Adoption of such a level gives greater certainty regarding the convergence of the solution.
When defining the stress-strain variables, one must specify the degradation variable d c , which is 0 for an undestroyed material and 1 for a total loss of the ability to carry stress: These values can also be read from the uniaxial compression tests by calculating appropriately the ratio of stress values of the declining branch of the graph to the compressive strength of concrete.Hence, the CDP model specifies strains as: In Eq. ( 5), E 0 denotes the initial modulus of elasticity of the undestroyed material.Based on the plastic strain values, it is possible to determine the stress σ c for uniaxial compression and their effective values  c describing the size of the flow and failure surfaces: The course of the material behaviour can also be determined without uniaxial compression testing, knowing only the average compressive strength of concrete f cm .Then, the second value necessary to start the analysis of the variability of the stress-strain curve is the modulus of longitudinal elasticity of concrete E cm .Its value can be estimated based on dependencies available in the code of practice [20]: Other values indicating the location of the characteristic points on the graph are the strains while reaching the strength ( c1 ) and the limit strains at failure ( cu1 ).
In turn, the tensile strength of concrete in the uniaxial tension state is rarely determined in a direct tension test due to implementation difficulties and large result variability.Therefore, one usually uses indirect methods, such as splitting samples or a beam bending method.In the absence of laboratory tests, it is possible to derive tensile strength from the compressive strength in accordance with [20]: The CDP model uses the notion of cracking strain  t ck .This allows to consider the phenomenon called tension stiffening effect [21].Concrete in tension is not treated as an elastic-brittle body, but one takes into account the effects of phenomena such as the aggregate interlock in a crack and the adhesion of concrete to steel in the segment between the cracks.This assumption has a reason during the analysis with a fuzzy cracking image.In this way, the decrease in the stress in concrete in the tension zone does not occur suddenly but gradually.The post-cracking strains are defined as the difference between the total strains and elastic strains of the undestroyed material: In turn, the plastic strains  t pl are calculated in the same way as for compression after defining the degradation parameter d t : In order to build the stress-strain curve of concrete in uniaxial tension, one must determine the form of the softening function (Fig. 5 [15]).[15].

Fig. 5. The stress-strain curves of concrete in uniaxial tension
In the absence of data, the stress can be reduced linearly down to zero from the moment of reaching the tensile strength in the case of total strains that are 10 times higher than at the time of reaching the value f ctm [15].A detailed description of this function, however, requires model calibration with experimental results in a specific analysis case [22]: In Eq. ( 13),  cr denotes the strain at concrete cracking.As the effect of stiffening can significantly affect the analysis results and additionally this relationship requires calibration during a simulation, it is recommended to use a modification of the formula by Wang T. and Tsu for the softening function: where n denotes a parameter that determines the softening rate.

Test specimens
The correctness of the adopted concrete model was verified on the basis of the results of experimental tests of continuous composite beams subjected to bending [23].
Test specimens were made from steel girders (European I-beam, IPN 360 type) connected with a concrete slab along its entire length using two rows of pin fasteners.The fasteners with the diameter of 16 mm and height of 75 mm were welded to the top flange of the beam with a 20 cm spacing.Such a connection was intended to ensure the rigidity of the joint in the entire load range.The primary reinforcement of the slab comprised twelve 10 mm smooth reinforcement bars, laid in two rows upside and downside, along the entire length, in groups of 6 bars with the spacing of 8 cm (the reinforcement degree was 2%).The transverse reinforcement consisted of stirrups made of 4.5 mm bars, arranged with a 20 cm spacing (Fig. 6).

Fig. 6. Tested beam parameters.
The places of application of external loads were chosen in such a way that the test piece included both the zone of negative moments causing tension of the upper fibres of the concrete slab and the zones of positive moments causing their compression (Fig. 7).The beams were loaded gradually by applying incremental concentrated forces in four full load-unload cycles.The first cycle was always carried out in the range from 0 kN until the appearance of the first crack.In the subsequent cycles, the force was increased by approx.200 kN until the maximum assumed load of 700 kN was reached.The maximum value was reached in each cycle by the gradual increase of the force, introducing intermediate load values.The load increase rate was on average 10 kN/min.Fig. 7. Test piece load scheme together with the slab damage measurement area.
All computational tasks were solved in Abaqus/Standard software using the incremental iterative Newton-Raphson method [15].The problems related to achieving the convergence of the solution, caused by the non-linearity of the material model, were solved by viscosity stabilisation μ [24].The size of the load increment was also reduced (0.01 -1E-12) and the maximum number of the load steps (max.12000) was increased.The selection of the parameter μ was made iteratively after analysing the size of its impact on the results of the task.Finally, μ=0.0001 was adopted, allowing one to solve the task in more than 1.200 load increments formed in approx.4.000 iterations.The analysis of the behaviour of the test configuration shows that such a value of the viscosity parameter makes it possible to reach a compromise between the computational magnitude of the task and the accuracy of the obtained results.

Geometrical models
Due to the nature of the structure being modelled, which consists of two different materials with clearly different geometries, the numerical model included different types of finite elements that best described the individual components of the beam.
Hence, the concrete slab was modelled with an 8-node brick element with reduced integration (C3D8R), and the I-beam with 4-node shell elements with reduced integration (S4R) and the primary reinforcement bars and stirrups with 2-node linear beam elements (B31).
Joining of the concrete slab with the upper flange of the steel beam varied depending on the distribution of bending moments along the length of the beam (Fig. 7).In particular, the joint was modelled in the tension zone as a discrete connection representing the occurrence of rigid connectors in the elements subjected to the laboratory tests.On the other hand, in the positive moment zones where the slab was compressed, the full joining of the upper surface of the steel beam with the lower surface of the reinforced concrete slab by means of a "continuous" (tie) connection was considered.

Identification of constitutive parameters of the concrete model
While identifying the constitutive parameters, it was assumed that concrete behaved in a linear-elastic manner up to 0.4•f c in compression and up to f t in tension.For compression, the elastic limit is the point described by the stress value of 28.06 MPa and the strain value of 0.0006450.Concrete hardens in compression from the stress value of 28.06 MPa to the value of 70.15 MPa.In this range, any strain relief occurs along a straight line, whose slope is E cm .Upon reaching stress that is equal to 70.15 MPa, concrete in compression softens, which is accompanied by the degradation of stiffness as well.After reaching the critical stress for concrete, the stiffness is reduced.
The stress-strain dependence of the non-linear behaviour of concrete under uniaxial compression (Fig. 8) was assumed in accordance with the recommendations [20]: The strain values corresponding to the greatest stress and the limit strain were determined from: Then, the variables were transformed into inelastic strains  c in , according to Eqs. ( 8)-( 9) and degradation variables d c were calculated as the ratio of stress values in the descending branch of the graph to the compressive strength of concrete.In this way, based on Eq. ( 12), plastic strains  c pl are determined during the numerical analyses.Similar deliberations were made for tension.In accordance with the earlier assumptions, it was assumed that concrete behaved elastically up to the stress value of 3.82 MPa.According to Eqs. ( 17)-( 18), cracking strains  t ck were determined, on the basis of which plastic strains  t pl are calculated taking into account the degradation parameter d t (the maximum concrete degradation equal to 0.99 was adopted).Based on the model calibration, the form of the softening function was determined with the assumption that stress can be reduced linearly from the moment of reaching the tensile strength (Fig. 9).
Table 1 shows the material parameters that determine the behaviour of concrete in a complex stress state.

Results and discussion
Identification of the slab cracking was done on the basis of the analysis of the maps of damage defined by the changes in the size of the DAMAGE parameter, i.e. the degradation of stiffness, which illustrates the material damage.
It should be remembered that the Concrete Damage Plasticity model of the material does not enable the formation of cracks in a discrete way taking into account the material spalling (its losses).It results only in the gradual exclusion of finite elements from the interaction.In this way, their "bonding" occurs, as well as their further involvement in the transfer of strain to the adjacent elements.This is one of the imperfections in the model, which has no significant effect on the behaviour of the entire test configuration, though.
The maps of damage defined by parameter d t in the tension zone and parameter d c in the compression zone, presented in the Fig. 10 and Fig. 11, can be identified with the spots where cracks appeared in the concrete slabs of the tested composite beams.It can be observed that the image of damage obtained in the numerical analysis corresponds qualitatively to the picture of the distribution of cracks obtained during the experimental tests.
The analysis of the maps of damage defined by parameter d t also allows for tracing the process of the formation and development of cracks with an increasing load.In the case of the analysed beams, the first damage to the slab concrete appeared at the support axis, with the load causing tensile stress in concrete that was equal to the tensile strength.In the initial phase, as the load continued to be increased, more cracks began to appear simultaneously on both sides of the support, spreading towards the centre of the span.The subsequent increase in the load caused the density of the damage zones (cracks) to increase.A similar phenomenon was observed during the empirical research [23].In addition, the verification of the numerical model was carried out based on the comparison of the obtained displacements of the tested beams and the numerical model (Fig. 12 and Fig. 13).The deflections were measured at two points where the external load was applied.In addition, readings of vertical displacements were taken in the longitudinal beam axis, above each support.
The analysis of the relationship between the load and the displacement confirms the good agreement of the solutions.The mean percentage differences between the computer calculation results and the measured values oscillate between ~5% and ~8%.
The graphs also show that the beam stiffness in the experiment and the stiffness of the numerical model are similar to each other.This confirms the correctness of the selection of degradation variables d (parameter d t , in particular) in the concrete model.concrete along with its progressive degradation as the stress values increase, however, can cause some problems related to achieving the convergence of the solution [24].These problems are mainly related to the tensile response.FEM techniques based on reducing the load increment size or increasing the maximum number of steps while solving the problem using the Newton-Raphson approach may prove insufficient.Hence, the Concrete Damage Plasticity model includes the viscosity parameter μ, which is used for constitutive equation regularisation (allows to slightly exceed the plastic potential surface in some sufficiently small steps).The very idea of viscoplastic regulation consists in such a parameter selection so that the ratio of the time step of the task to the value of μ would approach infinity.Such a regulation method forces the parameter μ to be selected multiple times to check its impact on the results of the problem and thus, determine its minimum value.
The results of the studies make it possible to formulate the following conclusions: -The material damage plasticity model is one of many possibilities for modelling concrete [25], which can be used for computer simulation of the behaviour in advanced load states of the entire structure and its individual elements.However, it is intended primarily for the analysis of reinforced concrete structures and composite structures of various types, such as steelconcrete and even concrete-concrete.
-Before the start of the analysis of a structure, one should determine the parameters of individual materials and constitutive relations between the strain and the stress.
-Due to the nature of the damage of concrete, some quantities should be adopted rather than determined as part of laboratory tests.Hence, the correctness of these assumptions should be checked by comparing other parameters, e.g. the deflection of the modelled structural element.Therefore, the results of the numerical analysis often require that the calibration of the model parameters be performed several times.
The author would like to thank the Gotowski Budownictwo Komunikacyjne i Przemysłowe Sp. z o.o.company for providing the testing materials free of charge.Laboratory studies were co-funded under the statutory measures of the Ministry of Science and Higher Education No. S-50/B/2012.

Fig. 2 .
Fig. 2. Deviatoric cross section of failure surface in concrete damage plasticity model and hyperbolic surface of plastic potential in meridional plane [15].

Fig. 10 .
Fig. 10.The comparison of crack patterns for bending beam in support section: a) beam with CDP numerical model, b) the fracture path, which is observed in experiment.

Fig. 11 .
Fig. 11.The comparison of crack patterns for bending beam in span section: a) beam with CDP numerical model, b) the fracture path, which is observed in experiment.

Fig. 12 .
Fig. 12.Comparison of the vertical beam displacements measured and calculated in support section.

Fig. 13 .
Fig. 13.Comparison of the vertical beam displacements measured and calculated in span section.