Model for the Simulation of the Time-Dependent State in RC Elements

: The paper presents an upgrade of the previously developed model for nonlinear 3D analysis of concrete structures extended with the possibility of simulation of the long-term effects (shrinkage and creep) under long-term static load. The origin model is based on the so-called multi-surface principle with modified Rankin criterion for dominant tensile influences (appearance and development of cracks) and the modified Mohr-Coulomb criterion for dominant compressive states (yielding and cracking of concrete). The material behaviour is described with elementary material parameters (modulus of elasticity, Poisson’s coefficient and uniaxial compressive and tensile strength of concrete) by standard tests. Sufficient accuracy along with a simple and effective description of the very complex behaviour of reinforced concrete structures, make this model advan-tageous. Creep and shrinkage are based on the procedure given by the fib Model Code 2010 and extended with a special extension for non-linear creeping. Two simple examples show the capabil-ities of the model, while a good agreement between numerical and experimental results indicates that the developed model can well describe long-term effects in reinforced concrete structures, and that the model is appropriate for standard engineering practice.


Introduction
When analysing the behaviour of concrete, it is very important to consider its timedependent behaviour. Concrete changes the stress and deformation states over time, i.e., it adapts to the external load. Experimental tests of behaviour under long-term loads require significant financial resources and long-term occupation of laboratories, measuring equipment and personnel. On the other hand, numerical models cost significantly less and can be used numerous times with different material characteristics, different loads and other parameters. But, of course, the used model has to be reliable and must produce reliable results.
When we speak about the 'time-dependent behaviour of concrete', we often think of creep and shrinkage, which are two main time-dependent physical properties of concrete.
Concrete shrinkage is the volumetric change of concrete structures caused by moisture loss due to evaporation. It is a time-dependent deformation, which occurs without the influence of external forces and causes a reduction of concrete volume.
On the other hand, creep can be defined as the long-term deformation of concrete under a permanent load. Generally, when a load is applied, pressure changes the shape of the concrete structure, and the strain occurs immediately along the direction of the applied load. But, during this time, deformation increases without the increase in load. We can say that concrete adapts to load or concrete creeps.
When the continuous load is removed, the part of the strain decreases immediately and another part over a long period (reversible creeping). In addition, some strain remains permanently on the structure. This physical phenomenon occurs at room temperature and is caused by several reasons: number and position of internal voids, rate of viscous flow of the cement-water paste, and, generally considered as the most important, seepage of water out of the cement gel due to external load and air humidity. The rate and magnitude of creep for most concrete structures are closely related to the drying rate, but it may also occur in massive structures where drying of the concrete is very small. In these structures, most of the creep is caused by external pressure, which forces the flow of the water absorbed from the gel. Creeping in concrete can result in a considerable increase of strains and usually needs to be considered during the design and analysis of reinforced concrete structures.
At low-stress levels, the creep of concrete can be described accurately by the linear viscoelasticity model, the so-called linear creep model. It has been experimentally determined that the linear functional relationship (linear creep) is only valid for the small-stress levels of concrete. When the relationship between stress of long-term static stress and concrete uniaxial pressure strength is greater than approximately 0.4, the linear relationship between the short-term and long-term deformations is lost and the progressive nonlinear increase deformation in concrete occurs, which is called non-linear creep [1][2][3][4][5][6][7][8].
When the stress of concrete with a long static load exceeds 80% of the uniaxial compressive strength of the concrete, creep often leads to the breakdown of concrete in time [9][10][11][12].
None of the offered models provides a complete description of the creep phenomenon. Each of the models has advantages, disadvantages and specifics of the application. According to Naaman [7], creep is defined as a time-dependent deformation higher than the initial deformation in concrete exposed to permanent load, while ACI Committee 209 [8] defines creep as a time-dependent increase in deformation in hardened concrete exposed to long-term loads. With the model proposed by Bažant-Baveja [9], the accuracy of estimation of deformation due to shrinkage and creep of concrete was increased compared to the ACI 209 model [8].
Of course, all the improvements usually reduce the simplicity of the model's mathematical definition and increase the number of input parameters. The user needs to decide the level of accuracy he needs, depending on the sensitivity of the problem being analysed. Typically, models are defined for a uniaxial stress state and then mapped to a multiaxial stress state. In these mappings and defining the functional dependence of uniaxial and multiaxial stress states, experimental studies of concrete creep have shown that the corresponding Poisson's creep coefficient at exposure to multiaxial compressive stress is not constant over time [10][11][12]. It depends on the multiaxial character of the stress state and is therefore not isotropic, and its initial value is lower than in the elastic region of the material. According to the literature [13], resulting from a twelve-year study of concrete creep exposed to uniaxial and biaxial stresses, the variation of the Poisson coefficient is negligible. It ranges from 0.32 to 0.27 for autogenous creep, while for total creep it remains at boundaries between 0.30 and 0.35, which justifies the mathematical formulations of creep and contraction defined with a constant Poisson coefficient.
Concrete creep has the characteristics of viscous deformations and occurs in both compression and tension zones. Although the tensile zone creep is less pronounced than in the compression zone, due to the lower tensile strength compared to the compression strength, creep in the tension is interesting in many practical situations [14], for example in assessing the possibility of opening cracks due to changes in temperature [15]. The studies reported in the literature [16] analyse the influence of the water-cement ratio on the creep development when the concrete is subjected to dominantly tensile stress. The development of creep deformation in reinforced concrete structures over time can lead to redistribution of stresses in the structure itself, whose neglecting, in combination with other adverse effects, could significantly impair the load-bearing capacity and stability of the observed structures.
Concrete shrinkage represents the total time-dependent decrease of volume caused by the change in humidity of the concrete mix. When deformations of concrete are prevented, shrinkage leads to the occurrence of additional stresses in the concrete mainly transferred to the reinforcement, which contributes to the reduction of shrinkage deformation [10,11,[16][17][18]. Shrinkage also reduces the development of continuous time-dependent displacement [19], thereby reducing the creep effect.
The main goal of this research is to develop and validate an efficient numerical algorithm/model and identify its parameters by fitting experimental data. The final objective is to create software for simulating full spatial concrete structures under long-term load. This software combines two models. The first is according to Ref. [8], which provides the prediction of the creep coefficient at any time, and the second is a 3D model for describing reinforced concrete structures that provide a stress-strain state at any time under a long time loading.

Introduction
The model used as the origin is the full 3D spatial model based on spatial brick elements. The model can describe the full spatial (triaxial) behaviour of concrete, with all dominant nonlinear phenomena important for a realistic simulation of concrete, such as yielding of concrete during compression and hardening and softening after yielding, as well as the occurrence and propagation of cracks in tension. A detailed description of the model, the equations, the functions and the parameters is presented in Refs. [20,21], so only a brief review of this model is given in this paper.

Numerical Model for Concrete for Dominant Tensile Influences
One of the most commonly used models for describing behaviour of concrete structures under dominant tensile stress is Rankin's material law. According to Rankin's law, a crack occurs when at least one of the principal stresses exceeds the tensile strength of concrete. In areas where other stresses are compressive, the tensile strength of concrete must be reduced ( Figure 1). In Figure 1, is the reduced tensile strength, ( = 1, 2, 3) are the principal tensile stresses in considered directions and is the uniaxial compressive strength of concrete. Figure 2 describes the mechanism of opening and closing of cracks in a concrete body (sampling points). In the beginning, concrete is considered an isotropic material with an isotropic stress-strain state (Figure 2a). After the first crack occurs, concrete becomes an orthotropic material whose axis is directed perpendicular to the crack plane i.e., towards the maximal tensile stress (Figure 2b).  After the first crack occurs, its direction stays fixed till the end of the simulation-for all of the following load increments. Subsequent cracks can appear only perpendicularly to the first one, and tensile softening is described by linearly reducing tensile stresses. The adopted model for cracked concrete follows the linear constitutive law with the adjusted modulus of elasticity [20]. The occurrence of the first crack softens the material, and the appropriate coefficients of the material matrix have to be reduced. The occurrence of the second and the third crack at the same integration point (perpendicularly to the first crack, i.e., perpendicular to the first and the second crack) further reduces the appropriate coefficients.
The fixed orthogonal crack model assumes that a crack can occur in a plane perpendicular to the principal tensile stress when this stress is higher than the tensile strength or the reduced tensile strength of concrete, depending on the area where the stresses have been calculated (Figure 1).
In the calculating procedure, if the value of the strain in each direction exceeds the prescribed maximum strain value in a Gauss point, the modulus of elasticity for that direction becomes 0, and material failure occurs in that direction. All mathematical formulations necessary for describing these procedures are defined in detail in [20][21][22].

Numerical Model for Concrete for Dominant Compressive Influences
The modified Mohr-Coulomb law (the elasto-plastic material model) is applied to describe the nonlinear behaviour of concrete in the dominant compressive stresses area. The so-called multi-surface model presentation [20] (Figure 3a), where the yielding surface is composed of six planes in the principal stresses coordinate system, is used. The graphical representation of the elastic, hardening and softening behaviour of concrete, defined as a function of total plastic strains, is shown in Figure 3b.

Numerical Model for Reinforcement Bars for Short-Term Load Condition
In the developed model, the reinforcement bars in the global concrete mesh are described by a 1D elasto-viscoplastic model, and their position is arbitrary in space ( Figure  4). These reinforcement bars are presented by a second-order space function determined by its projections [20,22]. The assumption is that the reinforcement transmits only the longitudinal tensile and compressive forces, so its behaviour can be approximated by a uniaxial stress state. In this model, a uniaxial elasto-viscoplastic model that includes potential hardening and elastic unloading is used [20][21][22].
In this work, the original model is upgraded with a simulation of creep and shrinkage, so now it is possible to simulate long-term loads on concrete elements/structures.

Linear Creep of Concrete
Creep can be defined as the progressive accumulation of the plastic strain in the concrete element under constant stresses (external load) over a long period of time. It occurs under a long-term (permanent) load when the element is exposed to stresses at high levels but below the ultimate material strength. The deformation rate corresponds to the material properties, time and temperature and the level of load applied to the structure.
Total strain at time t, in basic form, can be determined as: where: ( ) -total strain at some specific time , -average immediate mechanical (elastic) strain recorded at the moment after loading, ( , ) -creep strain at some specific time t, for element loaded at time , ( , ) -total strain caused by shrinkage at time t, determined on an unloaded specimen that begun to dry at time . The strains in a concrete member over time are illustrated in Figure 5. Total creep strain at a specific time is the sum of the basic and the drying creep components. For the sealed concrete without any moisture exchange with the ambient medium, only the basic creep component exists. If the moisture exchange with the surrounding medium is possible, drying creep also occurs.
For the calculation of the uniaxial creep deformation, the method of Glanville and Dischinger, based on the assumption that the creep rate function is in correlation with the current concrete uniaxial stress and time t, was used. In the most general form, this can be expressed as: For practical purposes, time is often divided into discrete time intervals Δ , with tn = t and = + Δ , so the incremental version of Equation (2) becomes: where: Some of the most used internationally recommended types of models for the prediction of creep strains without a need for specific laboratory tests are fib Model Code 2010 (2010) [23], ACI 209R-92 (1992) [8], AASHTO LRFD 2010, AS 3600, EC-2, among many others. These models are all empirical, and they vary widely in their techniques. As input, they all require main material properties and concrete age at the time of loading as well as certain intrinsic and/or extrinsic variables, for example, ambient moisture and mixture proportion.
In this paper, the fib Model Code 2010 was used [23]. The analytical form of the creep function is given with the following expression: with parameters [23]: ( , ) -part of the creep coefficient that represents basic creep, ( , ) -part of the creep coefficient that represents drying creep, ( ) -concrete elasticity modulus at load time, -concrete elasticity modulus after 28 days. The total creep coefficient is equal to the sum of its components: The creep coefficient components that appear in expressions (4) and (5) are determined according to expressions (6) and (7), [23]: In expressions (6) and (7) is the mean compressive strength of concrete in MPa after 28 days, coefficient RH is the relative humidity in %, is the perimeter of the member in contact with the atmosphere [mm] and is the area of the cross-section [mm 2 ]. The variable is the adjusted concrete age at the time of loading that can be calculated according to [23] as: with -temperature adjusted concrete age that can be calculated from: In the above expressions, Δ is the number of days where temperature (in ℃) prevails, while the constant in (8)  Since the spatial discretisation deals with the full spatial problem, the problem with the full tensor of deformations , , , , , is calculated with equal creep coefficients for all the deformation components. In this case, the creep coefficients, or the creep coefficients increments for individual deformation components, are calculated according to the previously shown procedure for a one-dimensional problem.
Finally, if the stresses in concrete do not change significantly, the deformations can be determined using the effective modulus of elasticity:

Nonlinear Creep of Concrete
It has been experimentally verified that the linear correlation between the short-term mechanical deformation of concrete and the long-term deformation (creep) of concrete is valid only when the stress levels in concrete are low. According to many authors, see for example [1][2][3], when the ratio between the stresses in concrete under a long-term load and the nominal pressure strength of concrete exceeds 0.40, the relationship between the shortterm and the long-term deformations becomes nonlinear and a progressive increase in the deformation, and the cracking of concrete occur ( Figure 6). Moreover, when stresses in concrete under a long-term load exceed 80% of the nominal concrete compressive strength, creep causes the collapse of the concrete element. There are various attempts at analytical descriptions of the nonlinear creep of concrete. Some models of nonlinear creep of concrete, based on different rheological models and/or on uniaxial experimental tests, can be found in [3][4][5][6].
The linear form of the creep function of Equation (5) In this paper, the empirical expression for nonlinear concrete cracking proposed in [20,24] was used. Considering the need for a simple practical nonlinear creep model, a simple function for the factor which depends on the current stress in concrete due to the and the compressive strength of concrete , has been adopted according to the slightly modified expression from [24]: As it can be seen in Figure 7, this model correlates well with the models given in [1,2,23,25]. This model is very simple and can be used up to the stress levels of 0.8 . Above this stress level, the model is not entirely reliable. In Figure 8, for 0.4 ≤ | | ≤ 0.6 , it is visible that the proposed expression (12)

Shrinkage
In this work, similarly, as for creep, the model described in the fib Model Code 2010, which predicts the main time-dependent shrinkage behaviour of plain structural concrete, was used [23]. This model is very good for standard engineering practices because it is based on the concrete compressive strength as a major parameter for estimating long-term concrete deformation properties. In reality, shrinkage and creep strains highly depend on parameters related to the microstructure and the composition of concrete, for example, the water/cement ratio, the properties of the aggregate, the degree of hydration, etc. Since engineers in practice often have only the compressive strength of concrete as the concrete parameter, it serves them as an indirect indication of the microstructural properties of concrete.

Numerical Procedure for the Time-Dependent Analysis
The model for the stress-strain analysis of reinforced concrete elements and structures under a short-term load, briefly described in Section 2 of this paper and more detail in [20][21][22], has been expanded by the simulation procedure for structures under a longterm load, considering the long-term effects in concrete: creep and shrinkage.

Numerical Procedure for Concrete Creep
Once the initial state (short-term analysis) of the structure is calculated, the long-term structural analysis begins. The time loop is set up, and for every time step, the current creep coefficient and the current shrinkage deformation are calculated depending on the stress level in the structural element.
Depending on the calculated stress level in each (Gauss) point, the material stressstrain relationship (i.e., the secant elasticity modulus of the material) is modified. The modification is performed according to the current creep coefficient that causes local concrete orthotropy in the area of some Gauss points. The new orthotropic matrix of the material which forms the stiffness matrix of the system is determined. Then, the structure is analysed with the new stiffness matrix, and the new state of stress and deformation is obtained as the result of the time-dependent behaviour of concrete. A correction of the material matrix is performed in the Gauss points and directions with both compressive stresses and tensile stresses that are smaller than the tensile strength of concrete, i.e., in those Gauss points and directions where cracks do not occur.

The Spatial State of Stresses and Strains
In the first step, when stresses, for a given combination of external forces, are within the initial flow surface, the following well-known relationship between strains and stresses applies: where is the effective modulus of elasticity modified with the Poisson's coefficient ( ): and is the initial modulus of elasticity and the modulus of shear, which, according to the linear theory, is calculated by: As the modulus of elasticity is corrected for principal stress directions in each Gauss point, it is necessary to introduce the orthotropic material model. The orthotropic directions coincide with the principal stress directions and become the major material directions.
Based on expression (14) for orthotropic materials in the spatial stress state, the relationship can be written as follows [27]: where: , and -moduli of elasticity in principal directions, , , , , and -Poisson's coefficients in planes indicated with numbers in indexes, , and -shear moduli in planes indicated by numbers in indexes Or shorter: The Poisson's coefficients can be calculated from the following expressions [27]: With the following condition, when ≠ [27]: = The shear modulus for the orthotropic material is calculated according to the following expression [27]:

The Correction of the Elasticity Modulus According to Creep
The numerical model performs an iterative calculation of the problem by time increments and load phases. At the end of each phase/time increment, stresses and strains in all of the Gauss points of the structure are known so that the principal strains and stresses (obtained from global stresses) and their directions can be easily calculated.
Based on the stress/strain state, the appropriate procedure is applied depending on whether the concrete is in the elastic region, the plastic (yielding) area, the region of tensile softening, or if the cracks have appeared in certain areas.
When the principal stress in a certain direction is known and the corresponding current deformation is obtained from the constitutive diagram of concrete, the current secant modulus of elasticity of concrete in that direction can be determined. For established directions of principal stresses ( , and ) , and their corresponding strains ( , and ), it is now possible to determine the secant current modulus of elasticity in the direction of principal stresses for each Gauss point with the following expressions: By determining the current effective (corrected) modulus according to expressions (10) and (22), their corrections are made concerning the development of the creep coefficient.
Based on expressions (10) and (22), the corrected modulus of elasticity for a single main direction and a single Gauss point with the compressive stresses and the tensile stresses smaller than the tensile strength of concrete at any given time can now be calculated.
This way, the corrected modulus of elasticity of concrete enters the matrix of the material, which now, for the spatial state of stress and a certain moment in time, takes the following form: Furthermore, the matrix of the orthotropic material should be transformed from the material coordinate system into the global coordinate system using the transformation matrix [20]. The orthotropic matrix of the material for the global coordinate system enters the finite element stiffness matrix. Now the creep deformations in principal directions can be calculated as follows: The current moduli of elasticity, as functions of principal stresses and strains, for one main direction are shown in Figure 9.

Numerical Procedure for Concrete Shrinkage
The numerical procedure for the simulation of concrete shrinkage is explained in detail in [26], so here only the necessary guidelines will be given.
The concrete shrinkage and temperature loads in concrete have very similar effects, causing deflections and strains but not stresses in elements that can be freely deformed [28,29]. In elements with constrained deformations only stresses occur, not strains. This fact was used for the development of the shrinkage numerical model based on the thermal analogy.
The first step in the numerical procedure is to determine the current stress-strain state of the element. After the time loop has been established, the total shrinkage for the observed time can be calculated from Equation (13). The shrinkage strain is a scalar value equal in all points for directions of all principal stresses. The second step is to calculate the temperature change T which formally replaces the total shrinkage strain: where ( − ) is the shrinkage strain, calculated from Equation (13), in the observed time and is the coefficient of the thermal expansion of concrete. This temperature change (25) in unreinforced concrete members causes only deformation , but not stresses (Figure 10c). If concrete members have reinforcement, part of the deformation is prevented , , so stresses occur (Figure 10d).
For this strain state, the appropriate stress state can be calculated: It should be emphasised that the stresses calculated according to expression (27) do not exist in reality and only serve to calculate the equivalent nodal forces from shrinkage (temperature). These nodal forces are returned to the global system, and the deformation state is calculated for them.
The newly calculated strain vector {ε} consists of the part representing the shrinkage {εs,sh} and the part caused by other external forces. The actual deformation state can be calculated by subtracting the previously calculated strains {εsh} from the newly calculated ones (Figure 10d) where vector represents the corrected strains, i.e., the real strain state of the structure. After the total strain vector has been obtained, the stress vector can also be calculated. To obtain the real value of the stresses, it is necessary to subtract the free shrinkage part of the stresses from the calculated value.
The corrected strains {ε'} and stresses {σ'} of concrete have to be used to control the occurrence and propagation of cracks in tension and yielding of concrete in compression [26,29].
The flowchart of the proposed numerical algorithm, previously described, is shown in the Figure 11.

Verification of the Numerical Model-Examples
This chapter offers three examples that served for the verification numerical model presented, i.e., the computer program Precon3D.

Three Concrete Walls-Long-Term Testing
The first example taken from the reference [31,32] was chosen. In these papers, the authors showed the results of seven tests and measurement of concrete shrinkage for three walls (ST1, ST2 and ST3), as shown in Figure 12. In the cited literature, the experimental results were compared with the results of the most important and most commonly used models for predicting concrete shrinkage: Model Code 2010, Eurocode 2, Model ACI 209-R92, Model B4 and Model B4s, which served to compare the results of all models with the model developed by the computer program Precon3D. The literature describes the program of experimental testing walls, from which all the necessary material parameters and geometry of the walls were taken to simulate the above experimental test with a computer program Precon3D. An experimental study focused on shrinking concrete by drying began in 2012 at the CTU Faculty of Civil Engineering in Prague. The test was performed on three walls of different thicknesses (b), 200 mm, 400 mm and 800 mm, while the height and width of the wall were the same and amounted to 800 mm ( Figure 12). Four sides of the wall are fixed (constraint of displacement), thus simulating a one-dimensional shrinkage of the concrete. Measurement of the strains in concrete began three hours after concreting at four measurement sites in the walls as follows: for two walls (ST1and ST3) for 522 days and for the third ST3, 300 days. Due to most common use, C30/37 concrete was employed, and the basic parameters of the concrete mixture were: water-cement factor (w/c) 0.50, cement content 360 kg/m 3 (Portland cement 42.5 R), aggregate content 1760 kg/m 3 and maximum aggregate size 16 mm. The compressive strength of concrete, measured on a specimen at the age of 28 days, was 55 MPa. Ambient conditions were carefully recorded in the laboratory every five minutes, and the measured mean relative humidity was 39.1%, while the average temperature was 21.5 °C for the period of 552 days. In the computer program Precon3D, walls ST1 and ST2 are discretised by three-dimensional (3D) 27-node Lagrange finite elements of 100 × 100 × 100 mm, while wall ST3 is discretised by the same finite elements of 200 × 200 × 200 mm, as shown in Figure 13. The strains of such modelled walls with given material parameters and early defined boundary conditions were monitored during 522 days for the two of the walls and 300 days for the third wall. The results obtained by the developed computer program Pre-con3D were compared with those from the literature obtained by experimental testing.
According to the results represented in Figure 14, it can be concluded that the numerical model Precon3D gives a very good agreement of the results with the experimental one. The biggest difference in results is obtained in the first 40 days of shrinkage, which can be attributed to the fact that the shrinkage model Fib Model Code 2010 gives a faster increase in deformation in the first days and then decreases over time. In the mentioned literature [31,32], the experimental results were compared with the results of the commonly used models for predicting concrete shrinkage: Model Code 2010, Eurocode 2, Model ACI 209-R92, Model B4 and Model B4s. So, we used all these results for comparison with results obtained by numerical model Precon3D individually for each wall type ST1, ST2 and ST3.
From the above-presented results, it is evident that the numerical model Precon3D not only gives good agreement with the experimental ones ( Figure 14) but generally gives better results in predicting shrinkage deformations than other models presented in the literature [31,32]. The word 'generally' is emphasised here because model EC2 gives better forecasting results in some phases of shrinkage (approximately from 28th to 220th days of monitoring) but before and after that period shrinkage prediction is less accurate in comparison with results obtained by Precon 3D (Figures 15-17). However, it can be seen that in the later phase of monitoring (after 220 days), the developed numerical model Pre-con3D gives the best agreement of the results with the experimental one, which can be very important for computer simulation of structural monitoring.   The numerical model Precon3D using the Code 2010 model to calculate shrinkage deformation gives better results than those obtained by the Code in the cited literature. It can be seen that the developed 3D nonlinear material model for concrete represents a very good and originally developed numerical model for 3D analysis of nonlinear behaviour of concrete with the inclusion of long-term strains.

Simple Supported Concrete Slab
To verify the developed model, a numerical simulation of the concrete slab presented in [25] was made. A detailed description of the experiment is given in [25], and only the most important data are repeated.
The main geometrical characteristics of the tested slab are given in Figure 18. The reinforcement bars were made from steel with the declared strength of 560 MPa, the declared yield point of 500 MPa and modulus of elasticity Es = 200 GPa. The concrete was prepared with the classic Portland cement and the lime aggregate with the maximum grain size of 8 mm, without additives. The slab was concreted on a vibrating table.
The concrete uniaxial compressive strength, the tensile strength and the modulus of elasticity were tested with the standard procedure, and the results were sorted in Table 1. After t = 90 days from the day concrete was placed, the slab was put on supports, exposed to self-weight and a total load of P = 10 kN. The load P was uniformly distributed over the central area of the slab 0.9 m × 0.9 m and applied in five equal increments, 2.0 kN each (one row of concrete cubes- Figure 19). After the full load P = 10 kN had been applied, it stayed for 360 days. After 3, 7, 14, 28, 90, 180 and 360 days from applying the full load, the slab deflection, the concrete compressive strain and the reinforcement tensile strain were measured at selected points for each load increment. The picture of the fully loaded tested slab is shown in Figure 19.

Numerical Analysis by Precon3D and the Comparison of the Obtained Results
In Figure 20, the finite element mesh for the concrete slab is shown. The reinforcement is shown with thick red (x-direction) and thick blue (y-direction) lines. The properties of each material are chosen according to values from the experiment. Some results from the numerical model are shown in Figure 21. Figure 21a-c show the measured slab deflections during loading (90 days after concreting) and the comparison with the numerical results is made. In Figure 21a,b, the comparison of the numerical and the experimental results for measuring points MU-1 and MU-2 is shown and a good agreement of the results is visible. Figure 21c shows the comparison of the numerical and the experimental results for measuring point MU-3 and, unlike measuring points MU-1 and MU-2, it shows a slightly lower agreement of the results, especially in the region of the appearance of nonlinearity. At this measuring point, a lower load level causes the appearance of more significant nonlinearity than in measuring points MU-1 and MU-2. The most pronounced difference between the numerical and the experimental results are in this region. There are many reasons for this difference, but the most prominent reason is the high amount of reinforcement placed in the y-direction. This significantly prevents the development of shrinkage strain that causes the appearance of higher tensile stresses in concrete in the y-direction than in the x-direction. These tensile stresses cause an earlier appearance of cracks in concrete, i.e., nonlinearity. For the final load level, there is a good agreement between the numerical and the experimental results.
A significant increase in time-dependent deflections of the slab after the full load P is applied can be observed in Figure 21d-f for the experimental and the numerical results. This results from nonlinear creeping caused by high-stress levels in the slab and the further spreading of cracks in the concrete tensile zone. The numerical results follow the measured results very well for all measuring positions.

Reinforced Concrete Girder
The second example is the analysis of a reinforced concrete girder taken from Ref. [33]. According to this reference, this girder with all of the ambient influences was monitored for two years. The geometry, the position of the load, the measuring line and the material characteristics are presented in Figure 22 and Table 2, all of the data were taken from Ref [33].    Figure 23a shows the state of stresses resulting from the shrinkage of concrete before the application of the external load (day 29). From the figure, it can be seen that the tensile stresses appear in the zone where reinforcement is placed, while the compressive stresses appear in the zone without reinforcement. It also can be seen that the tensile stresses are at the level of the tensile strength of concrete and that the first crack was formed even earlier and at this point, new cracks are opening. The compressive stresses are much smaller than the compressive strength of concrete. Figure  23b shows the stress conditions after the external load is applied (day 30), and a significant increase in the compressive stresses is visible. On the other hand, the tensile stresses remain at the level of the tensile strength of concrete which indicates the opening of new cracks and the widening of old. The shown crack condition (c) corresponds to the maximum size of the external load (day 30). Here the compressive stresses have reached their maximum value of 13.42 MPa, while the cracking of concrete in the tensile zone is considerable. Thus, the concrete loses its load-bearing capacity under tension, resulting in tensile softening of the concrete and a decrease of the tensile stresses in concrete. Figure 23d shows the stress distributions for a permanent external load, i.e., the changes in stresses only due to the long-term effects: the creep and the shrinkage of concrete. It can be seen that both the tensile and the compressive stresses in concrete are reduced. The decrease in the compressive stresses is due to the lowering of the neutral axis towards the tensile zone caused by the shrinkage of concrete, while the decrease in the tensile stresses is due to the increase in the cracking of concrete in the tensile zone. This is caused by the redistribution of stresses in the girder as a consequence of long-term deformations of concrete. The reduction of the compressive stresses at the end of the observation (day 730), relative to the stresses at the time of the application of the full external load (day 30 Girder length(m) the reduction of the tensile stresses is about 56%. The dominant stresses in the x-axis are analysed in more detail, while the normal stresses in the y-axis ( ) are caused by the triaxial stress state and are negligible compared to the normal stresses in the x-axis, and their change over time is also negligible. The normal stresses in the z-axis ( ) are also a consequence of the triaxial stress state, and are also negligible relative to the stresses in the x-axis, although they are slightly larger than the stresses in the y-axis, which is clearly due to the external load in the z-axis. These stresses decrease over time due to the creep and the shrinkage of concrete. The compressive stresses at the end of the observation period (day 730) are about 14% lower than at the time of the application of the full external load (day 30), while the tensile stresses on day 730 are about 12% lower compared to day 30. Furthermore, the analysis of the stresses in the reinforcement bar at the same moments in time as for concrete is given. Figure 25a clearly shows that just before the external load is applied, only the compressive stresses, as a consequence of shrinkage, occur in the reinforcement. These stresses in the reinforcement cause the appearance of compressive forces in the reinforcement bar which are transferred to the concrete as tensile forces and cause the stress states in the concrete shown in previous Figures 23 and 24. In Figure 25b it can be seen that, immediately after the application of the external load, the stresses in the reinforcement bar change from compressive to tensile. This indicates that the concrete has cracked due to shrinkage and that after loading, all of the load-bearing capacity of the tensile zone was taken over by the reinforcement. Figure 25c shows the state of stresses in the reinforcement bar at the time of the application of the maximum size of the external load, and a considerable increase in stresses is seen compared to those that were present at a smaller load size. It points to the opening of new cracks and the expansion of the existing ones. In Figure 25d, it is visible that the propagation of cracks due to concrete creep will cause an increase in the stresses in the reinforcement bar. The maximum stress in the reinforcement bar on day 730 is about 28% higher than it was on day 30 after the application of the full external load.
This increase in tensile stresses in the reinforcement bar is shown in Figure 26, where the stress increase at the reinforcement bar point near the middle of the span of the girder, is visible between day 30, after the application of the full external load, and day 730. The figure shows that the increase in these stresses, caused by concrete creep, is greatest in the first 10 days of creeping, while in the later phase of creep it stabilises and the increase in stress is much smaller. So, in the later stages, due to creep, the stresses in the reinforcement bar reach a constant magnitude.

Conclusions
The paper presents a developed numerical model for the analysis of spatial reinforced concrete structures with emphasis on their behaviour under long-term loads. The previously developed model for short-term loads has been extended by the possibility of simulating long-term effects in concrete: creep and shrinkage. Thus, a modified nonlinear material law using the 3D representation of Rankine's material law for dominant tensile stress and a multi-surface presentation of modified Mohr-Coulomb material law for dominant compressive stresses has been extended with the phenomenological model proposed by fib Model Code 2010, with modification to include nonlinear creep. Analysis of the influence of shrinkage on reinforced concrete structures was performed using an equivalent temperature load. For the observed type of construction of appropriate material and geometric characteristics, using the adopted shrinkage model, the shrinkage deformations of concrete (without reinforcement) are calculated for the directions of main stresses. As a ratio of the deformations thus obtained and the coefficient of thermal conductivity of concrete, equivalent nodal temperature forces are obtained with which the actual system is then loaded (considering the reinforcement).
The analysis of the effect of creep on the behaviour of concrete uses the model of correction of the modulus of elasticity of concrete in the direction of the main stresses. The moment of application of external load represents the starting moment for the creep coefficient calculation. The creep coefficient depends on stress level, and it is calculated according to expression (12), for each direction independently at every time step and load increment step. The stress and strain states at each Gaussian point for the direction of the principal stresses are calculated for each moment, so the modulus of elasticity of concrete is known at each observed moment. If some of the principal stresses are lower than the tensile strength of concrete, the modulus of elasticity is corrected in the direction of that stress considering the magnitude of the creep coefficient at that moment. This calculation is performed for each Gaussian point so that by correcting the modulus of elasticity, orthotropy of the material at the level of each Gaussian point is achieved. With such a corrected modulus of elasticity, the concrete stiffness matrix is entered, and a new concrete stiffness matrix is obtained. By using the stiffness matrix thus obtained, the strains caused by creep are obtained. Relating to time analysis, the magnitude of the time step does not have to be constant. However, as the increment of creep and shrinkage deformations is initially higher and decreases with time, a shorter time step is taken at the beginning of the observation, which increases in the later phase of the calculation. This interpretation of creep and shrinkage through the numerical model Precon3D resulted in a model that can conduct a time-dependent analysis of the behaviour of reinforced concrete structures.
Verification of the presented model is shown through the analysis of three different examples: analysis of concrete walls, analysis of reinforced concrete slab and reinforced concrete beam. All structures were monitored and analysed for a long time with and without the influence of external load. The examples were selected to allow verification of the developed model for different structures and loads.
The analysis of the first example verified the accuracy of the calculation of long-term strains in concrete. A comparison was made with experimental measurements and the results obtained by calculating other most commonly used models: Model Code 2010, Eurocode 2, Model ACI 209-R92, Model B4 and Model B4s. It can be seen that, especially in the later phase of monitoring, the developed numerical model Precon3D gives the best agreement of the results with the experimental one, which can be very important for computer simulation of structural monitoring.
In the second example, a reinforced concrete slab was loaded 90 days after concreting was analysed, and its behaviour was monitored for the next 365 days. The analysis of this example with the Precon3D program shows a very accurate calculation of nonlinear creep caused by high-stress levels in the slab and further propagation of cracks in the tensile zone of concrete. Moreover, the numerical results obtained by Precon3D follow the measured results very well for all measuring positions.
The third example analyses the behaviour of the reinforced concrete beam, where the influence of shrinkage on the behaviour of concrete and reinforcement is analysed. The results obtained in the numerical simulation by Precon3D agree with the theoretical logic of the influence of shrinkage on the behaviour of concrete and reinforcement. It has been noticed that shrinkage has a certain influence on the behaviour of reinforcement, so it should not be neglected in the long-term analysis. As concluded from the analysis of numerical results, the concrete shrinkage deformation represents the nature of concrete behaviour and cannot be treated as 'classical' deformation (does not cause stress change in concrete), such as that caused by an external load.
The Precon3D numerical program is designed to leave users with the choice of conducting short-term or long-term analysis. If it is a long-term analysis, then there is a possibility to analyse only shrinkage, only creep, or both creep and shrinkage of concrete. In the analysis of concrete shrinkage, it is possible to perform the total shrinkage analysis as the sum of autogenous shrinkage and shrinkage during drying or only the analysis of shrinkage by drying.
Most computer programs/numerical models for 3D nonlinear analysis of reinforced concrete structures for short-term and especially long-term analysis require many specific parameters, so the main advantage of this model is that it can perform 3D nonlinear material analysis for short-term and long-term loading of reinforced concrete structures using basic material parameters obtained by uniaxial testing. The presented numerical model is reliable and applicable for monitoring structures even at high-stress levels, i.e., nonlinear creep.