A multiscale model for reinforced concrete with macroscopic variation of reinforcement slip

A single-scale model for reinforced concrete, comprising the plain concrete continuum, reinforcement bars and the bond between them, is used as a basis for deriving a two-scale model. The large-scale problem, representing the “effective” reinforced concrete solid, is enriched by an effective reinforcement slip variable. The subscale problem on a Representative Volume Element (RVE) is defined by Dirichlet boundary conditions. The response of the RVEs of different sizes was investigated by means of pull-out tests. The resulting two-scale formulation was used in an FE2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} analysis of a deep beam. Load–deflection relations, crack widths, and strain fields were compared to those obtained from a single-scale analysis. Incorporating the independent macroscopic reinforcement slip variable resulted in a more pronounced localisation of the effective strain field. This produced a more accurate estimation of the crack widths than the two-scale formulation neglecting the effective reinforcement slip variable.


Introduction
Durability of reinforced concrete structures is closely linked to the brittle nature of concrete itself, which makes them vulnerable to ingress of harmful substances, very often causing corrosion of reinforcement, cf. the works of Nilenius et al. [26][27][28]. Crack width is an important factor influencing the durability of the structure. Hence, crack widths are often limited in design codes [6]. To be able to model crack growth kinematics of steel rebars and their bond with concrete in an enriched solid finite element. Ibrahimbegovic et al. [16] represented the concrete strain by using embedded discontinuity (ED-FEM, cf. the work of Mosler and Meschke [24]), while describing the slip strain with XFEM and combined them both into a single solid finite element. A macroscopic beam element with higher order interpolation of the displacement field, able to describe crack opening and rebar slip was developed in [7]. Another one-dimensional finite element model of a reinforced concrete beam accounting for interaction between the cracked concrete and reinforcement, while keeping track of the slip between the constitutive materials and different length scales, was established by Floros and Ingason [13].
Due to its heterogeneous microstructure, concrete itself evinces complicated composite behaviour. Modelling each phase in detail in the whole structure would result in very large and computationally expensive models. One possibility to tackle this problem is to employ multiscale modelling methods for concrete. A range of multiscale methods suitable for modelling plain concrete has been extensively reviewed by Unger and Eckardt [35]. A few different multiscale models of plain concrete at mesoscale were studied in the literature. For example, Nilenius et al. [27] used mesoscale models to study diffusion phenomena in plain concrete. Wu and Wriggers [37] developed mesoscale models to study diffusion-thermal-mechanical coupling, while the damage behaviour in plain concrete was studied at mesoscale by Wriggers and Moftah [36]. Furthermore, a multiscale model for plain concrete linking macroscale, mesoscale, and microscale was developed by Nguyen et al. [25]. Within the range of multiscale modelling methods, Belytschko and Song [1] identified hierarchical and concurrent methods. The latter group is characterised by bidirectional link between the scales. Concurrent multiscale methods have been used to model reinforced concrete frames and beams by Sun and Li [33] and Sun et al. [34]. Such methods require fine resolution of damaged parts of the structure. Although viable for frames, for structures that display distributed cracking, e.g., in the presence of reinforcement grids, concurrent multiscale methods would require fine resolution of the substructure almost everywhere. On the other hand, hierarchical multiscale methods involve transfer of information from the fine scale to coarse scale. An example of such method is FE 2 approach, introduced by Feyel [9,10]. This method couples the scales in a nested way, as the response on the macroscale is obtained through computational homogenisation of the response of the subscale boundary value problem on a Representative Volume Element (RVE) located at the macroscopic quadrature points. Although also computationally demanding, this approach is very well suited for parallel computing, since all RVE problems are uncoupled and can, hence, be solved concurrently.
However, only a few works considered detailed bond-slip modelling with hierarchical multiscale methods. Lackner and Mang [18,19] have used multiscale methods in calibration of a macroscopic constitutive model for reinforced concrete. The model rested on the assumption that the fracture energy of concrete related to the opening of primary cracks should be increased to reflect the bond slip between steel and concrete. Analyses at reinforcement bar scale were used to calibrate the model parameters at structural member scale, and the model was applied in modelling of several real structures [19]. As a result, the failure mode and general behaviour of the reinforced concrete structure were well captured. In a recent work by the authors [32], a two-scale model of reinforced concrete was developed and used within FE 2 setting. Large RVEs were required to give acceptable results in terms of the crack widths, which was attributed to the model assumption of the reinforcement slip varying only locally, i.e., at the subscale. As shown in [16], ideally the whole rebar should be modelled in order to obtain accurate information about crack spacing and widths, due to slip and bond stress transfer along the bar. Modelling the whole reinforcement bar at macroscale may not be feasible in FE 2 . Instead, the reinforcement slip could be considered as a macroscopic variable and thus its transfer between macroscopic elements could be activated.
This approach is applied in this paper, where the existing two-scale model of reinforced concrete [32] is extended to allow for reinforcement slip transfer between the macroscale elements. The multiscale formulation is devised using the variationally consistent homogenisation, proposed by Larsson et al. [21]. The relevant equilibrium equations are established for a subscale RVE, and first-order homogenisation is used to obtain the large-scale response. Both concrete and steel possess an independent macroscopic component. It should be noted that even though the formulation maintains generality in three dimensions, the multiscale model is in this work used in two dimensions.
The remainder of the paper is organised as follows: After some preliminaries concerning the fully-resolved problem in Sect. 2, the corresponding two-scale formulation is derived in Sect. 3. Section 4 presents the subscale response of a reinforced concrete RVE subjected to varying reinforcement slip. In Sect. 5, the developed model is applied in a two-scale analysis of a reinforced concrete deep beam. The paper is concluded with Sect. 6, which contains some final remarks and an outlook to future work.

Fully-resolved problem
In this section, a boundary value problem comprising a reinforced concrete structure is considered, and the variational format of governing equilibrium equations is presented. For a more detailed derivation, the reader is kindly referred to Fig. 1 A two-dimensional reinforced concrete structure. For each reinforcement bar, longitudinal and transverse unit vectors e l and e ⊥ are defined Sciegaj et. al. [32]. A symbolic representation of a twodimensional reinforced concrete structure is given in Fig. 1.
The problem domain Ω = Ω c ∪ Γ int comprises the concrete (Ω c ) and reinforcement (Γ int ) parts. It is assumed, that the reinforcement does not cross the external boundary Γ ext := Γ u ∪ Γ t . For the 2D problem, displacement fields pertinent to concrete and steel are denoted u c and u s , respectively. Along Γ int it is possible to split u s and u c into components that are parallel and perpendicular to Γ int , i.e., It is noteworthy that the directions e l and e ⊥ need not be constant throughout the whole structure, but may be a function of the position. Similarly, as the unit vectors e l and e ⊥ are associated with each bar (and can vary along it), the orthogonality of the reinforcement bars is not required. Taking t c to be the thickness of the structure, b is the body force and σ c is the stress in the concrete phase, the equilibrium can be stated in strong form as follows For steel, it is assumed that the reinforcement can sustain both normal force and bending moments, i.e., it can be subjected to both longitudinal and transverse loads. In the former case, we have the normal force N s linked to the bond stress t Γ (distributed around the circumference of the bar S s ). The latter mechanism couples the bending moment M s and the transverse distributed load λ. Therefore, for each reinforcing bar equilibrium can be expressed as Next, we consider the steel/concrete interface as presented in Fig. 2, where the reinforcement bar was fictitiously protruded out of the concrete. Summing all the forces acting on the concrete along Γ int , the equilibrium condition for the interface can be found to be Lastly, it is assumed that there is no relative motion between the steel and the concrete in the transverse direction (no normal displacement jump), i.e., we introduce the following interface constraint: By adopting this formulation we assume the contact deformations in the transverse direction to be negligible. To maintain generality, only implicit (algorithmic) definitions of the constitutive relations are considered. Hence, evolving internal variables are considered, but omitted in the abstract notation. Thus, for both the concrete and the steel we consider the implicit relations Fig. 2 Steel-concrete interface. Boundary forces (normal forces, shear forces and bending moments) on the rebar cut out, as well as the boundary tractions on the concrete are omitted for clarity where ε = [u c ⊗ ∇] sym denotes the strain and s = u s,l − u c,l is the reinforcement slip. In order to state the variational format, (3), (4), and (6) must be recast into weak forms. Multiplying each equation by a suitable test function, integrating over the domain (and employing (5) in (3) together with neglecting body forces), we arrive at the definition of the fully-resolved problem: Find for suitable trial sets U c , U s,l , U s,⊥ , L defined as and the test space U 0 c , defined as The coupling terms indicated in the system (11)-(14) are defined as: The following forms are introduced pertinent to (i) Concrete: (ii) Bar action of the rebars: (iii) Beam action of the rebars:

Two-scale problem
The classical terms macroscale and microscale are commonly used in multiscale modelling literature. In this paper, reinforced concrete is treated as a continuum material with distinct reinforcement bars. Even though this is considered to be the "finer" scale (for the purpose of the model), the physical length scale is clearly macroscopic. In order to make the formulation more general and to better emulate the nature of the problem, the classically used terms macroscale and microscale were substituted with the terms large-scale and subscale, respectively. Furthermore, the terms macroscale and large-scale are used in the following interchangeably.

Variationally consistent homogenisation
The scale separation technique used in this paper is referred to as variationally consistent homogenisation, and is treated in detail in Larsson et. al. [21]. Employing the Variational MultiScale (VMS) ansatz, the unknown fields u c , u s,l , u s,⊥ are separated into the smooth (large-scale) and fluctuating (subscale) parts: where the superscripts M and s denote large-scale and subscale components, respectively. A consequence of the mentioned separation is a necessity to solve two problems instead of one. Namely, we have a large-scale boundary value problem defined on globally "smooth" fields, and a subscale boundary value problem defined on the fluctuation fields within an RVE. Since in practice, the integrals are numerically evaluated at the integration points, it suffices to solve the subscale problem only there. The RVE is then defined in a region Ω in the vicinity of large-scale integration point located atx. It is assumed, that the reinforcement bars do not change direction inside the RVE, i.e., the unit vectors e l and e ⊥ are constant for each reinforcement bar in the RVE.
In the two-scale model setting, the local field is replaced by the homogenised field, i.e., at each locationx ∈ Ω the field is approximated by the volume average on Ω (x). More specifically, for given functions f Ω and f Γ defined on Ω c and Γ int , respectively, we have where the subscale average f is defined as: Utilising the subscale averages, we can express the fullyresolved problem defined in Eqs. (11)- (14) as where we introduced the RVE-forms a ,l (u s,l , u s,⊥ ; δu s,l ) : a ,b (u s,l , u s,⊥ ; δu s,⊥ ) : To obtain the variation of the large-scale "smooth" field within the subscale region Ω , Taylor expansion is used. In this work, the "smooth" field is expanded up to the linear term, i.e., first-order computational homogenisation is used, cf. the work of Geers et al. [15]. We introduce two independent fields-one describing the displacement,ū, and one representing the reinforcement slip,s. The prolongation conditions of the "smooth" large-scale fields in the subscale region can therefore be expressed as Note that due to the assumed constraint (6), there is no deformation between concrete and steel in the transverse direction. This in turn prevents the macroscopic slip fields from enriching the prolongation condition for u M s,⊥ . Even though the physical interpretation of the macroscopic displacementū is intuitive, the meaning of the macroscopic slip,s, is further elaborated on in "Appendix A".

Large-scale problem
Splitting the test functions in Eq. (31) according to the VMS ansatz, we divide the fully-resolved problem into large-scale and subscale problems. The former is described with Assuming that the local field is sufficiently smooth at the external boundary Γ ext , we have l c δu M c ≈ l c (δū), and utilising the prolongation conditions (37)-(39) we obtain The quantitiesσ ,τ b ,σ s are named the effective stress, effective transfer stress and effective reinforcement stress, respectively. Finally, the large-scale problem is defined as: with the suitable trial and test spaces In two-dimensional setting,σ represents the membrane stress in the reinforced concrete, i.e., its unit is force/length. The physical interpretation of the effective transfer and reinforcement stress is presented in "Appendix A". The details concerning definition and implementation of the large-scale finite element possessing two independent fields can be found in "Appendix B".
Remark By divergence theorem, we have alsō where the summation over the discrete forcesR L ,R ⊥ and R M is performed at locations where the reinforcement intersects Γ , cf. Fig. 3. These forces together with the boundary tractiont are defined as follows:

Fig. 3 Discrete forces and traction at the boundary Γ of an RVE
with the scalar e ln defined as

General formulation
Following the VMS ansatz, i.e., the fully-resolved problem (31) is tested with functions pertaining to the fluctuation fields. In general, before specifying the boundary conditions on any given RVE, the subscale equilibrium can be expressed in the weak format as with the test spaces defined as follows: Furthermore, in the general case (before boundary conditions have been specified), the boundary terms are given as where the discrete forcesR L ,R ⊥ ,R M , and the tractiont were defined in Fig. 3. Note that, although the above formulation maintains generality, the boundary conditions on the local fields need to be specified in order to produce a solvable system. Even though a couple of choices are possible, the focus of this paper is put on Dirichlet boundary conditions, as they have numerous advantages, such as reliability and ease of implementation [32]. For brevity, we define the macroscopic strain,ε, and macroscopic slip gradient,ḡ, as the gradients of the macroscopic displacementū and slips, respectively: Note that, in general, the gradient of the slip does not need to be symmetric.

Dirichlet-Dirichlet boundary conditions
One of the conspicuous choices is to prevent any fluctuation fields from forming at the boundary of the RVE, Γ . This way we restrict the total field to be equal to its macroscopic part, according to the VMS ansatz: Employing the Dirichlet boundary conditions in (55)-(58), the Dirichlet-Dirichlet (DD) subscale problem takes the for suitable trial sets and test spaces

Subscale response
The effective response of an RVE of reinforced concrete subjected to macroscopic strain was studied by the authors in [32]. The new formulation in the current work requires that, in addition to the macroscopic strain, also the macroscopic slip and slip gradient are imposed upon the RVE. In order to study the response of reinforced concrete within the developed framework, a series of pull-out tests was simulated. In these analyses, the slip of the reinforcement grid was gradually increased, while the concrete stayed intact, therefore resembling a conventional pull-out test [4] where a rein-

Modelling choices
For the RVE, a two-dimensional arrangement of a periodic reinforcement grid (φ 1 bars longitudinally and φ 2 bars transversely) was considered. To investigate the sensitivity of the developed two-scale framework to the size of the RVE, three sizes of unit cells were studied, all with a shape of square with side L , and thickness t, see Table 1. In all unit cells, spacing of the reinforcement bars in the grid was equal to 0.2 m. Type B500B [3] reinforcement was modelled with beam elements and Von Mises elastoplastic material model with strain hardening.
For the bulk of the RVE, bilinear quadrilateral elements were used to model concrete of grade C30 [12]. The Mazars model [22,23] is a commonly used isotropic continuum damage formulation for concrete. In this paper, the damage evolution law in tension, g t (κ), was modified so that the tensile stresses disappear after reaching a certain level of cumulative equivalent strain κ, i.e., where ε 0 = f ct /E c is the concrete strain at the onset of softening, and ε f is a mesh dependent parameter depending among all on the fracture energy, G F , and the element size, L e , cf. [29] for more details on the alteration of the damage law. After calibration of the model for quadrilateral elements of size L e = 0.02m, the parameters ε 0 = 8.63 × 10 −5 and ε f = 2.5 × 10 −3 were chosen. The model parameters pertinent to the damage evolution law in compression, A c = 2.6, B c = 800, are given for the sake of completeness.
Interface elements with thickness equal to the circumference of the rebar were used to model the steel/concrete interface according to the bond-slip law included in Model Code 2010 [12]. Furthermore, the interface constraint in Eq. (6) was satisfied by tying the translational degrees of freedom in the steel, to the corresponding degrees of freedom in the concrete along the bars. Material parameters for the concrete, steel and the interface are presented in Table 2, and the constitutive relations for all materials are schematically depicted in Fig. 4. Note, that in this example similar material properties for steel, concrete and the interface were chosen for all elements. However, they could all be varying, i.e., different properties could be present in different parts of the model.

Effective response at reinforcement pull-out
In the simulations, the macroscopic slip (in one direction) was prescribed on the RVEs in steps of Δs = [5 × 10 −5 m, 0] T , while the macroscopic strain in the concrete was kept constant atε = [0, 0, 0] T . Similarly, the macroscopic slip gradient was also kept at zero, i.e.,ḡ = [0, 0, 0, 0] T . After imposing the macroscopic quantities and solving the subscale problem, the effective transfer stress,τ b , was computed according to (48). From this homogenised value, the average bond stress per reinforcement bar in the RVE could easily be recovered, see "Appendix A". The effective relation between the imposed macroscopic slip,s, and the average bond stress (effective bond stress),t Γ , is reproduced for all RVEs in Fig. 5. From Fig. 5, it can be inferred that for the smallest unit cell (1 × 1), the bond-slip law for the interface was directly recovered, as both bond stresses and slip values corresponded to the input parameters. The same cannot be said about the larger RVEs, where a greater macroscopic slip was needed to be prescribed to extract the full bond-slip behaviour (up to the plateau). Furthermore, the average bond stress peak values for 2 × 2 and 3 × 3 RVEs were decreasing. The reason is that longer reinforcement bars did not have a constant distribution of bond stresses along the bar. In fact, Cairns and Plizzari [4] distinguish between "short" bond lengths (smaller than 6 times the bar diameter), for which the bond stresses can be considered constant along the bar, and "long" bond lengths Based on the results, it can be educed that the situation in the smallest RVE corresponded to the short bond length (although being 10 times the bar diameter; i.e., a bit longer than 6 times its diameter). Thus, the slip and bond stresses distributed almost uniformly along the bar. The largest bond stress was reached along the entire bar, and therefore the peak value of the effective bond stress was equal to the peak value of the bond-slip law. The rebars in the 2 × 2, and 3 × 3 unit cells were longer than 15 times their diameters, and thus pertained to long bond lengths. Neither slip nor bond stresses were constant along the bars. Although the maximum bond stress was reached in the middle part of the bar, the ends of it experienced larger slip values and consequently lower bond stresses. Therefore, the effective bond stress was in both cases lower than the peak value of the bond-slip law. The longer the bar, the lower the average (macroscopic) bond stress.
To summarise, the bond-slip law was recovered from the pull-out test either directly (for short rebars) or indirectly (for longer rebars). The discrepancy between the effective bondslip curves was caused by the nonuniform distribution of bond stresses along the bar, for varying embedment lengths. It is noteworthy, that prescribing the macroscopic slip via Dirichlet boundary conditions, i.e., prescribing values only at the boundary of the RVE, directly affects primarily the end part of the longer reinforcement bars. Both the amplitude and fluctuation of the slip are highly dependent on the size of the RVE (which translates directly to the length of the rebar). For a chosen size of the RVE, the slip distribution is more or less "resolved" in the subscale model. Therefore, the interpretation of the macroscopic slip fields depends on the RVE size. A practical implication of this is that subscale results must always be consulted in order to get the correct interpretation of the large-scale results. On the other hand, it could be argued that a better way of imposing the reinforce-

Structural response
For this numerical study, a simply supported deep beam of reinforced concrete was considered, according to Fig. 7. The same subscale features (thickness and reinforcement layout) as described in Sect. 4 are assumed uniformly distributed in the deep beam. Since the thickness of the structure (0.2 m) was much smaller than the other dimensions, a twodimensional model in plane stress was deemed sufficient for modelling the beam. The reinforcement comprised a uniform grid made of φ20 bars longitudinally and φ8 bars transversely, spaced 200 mm in both directions. The loading and support platens (0.76 and 0.9 m wide, respectively) were used for applying the load. The purpose of the study was to analyse the structural response of the deep beam, in terms of its load-deflection relation, crack widths and crack pattern. Furthermore, it was of interest to see how well the developed two-scale formulation with macroscopic variation of slip (Sect. 3) performs in comparison to a single-scale analysis in full-resolution, and to the two-scale formulation disregarding the reinforcement slip at macroscale, developed by the authors in [32]. For the computations, an open source C++ code OOFEM (www.oofem.org) [30] was used.

Single-scale model
The fully-resolved model comprised the concrete (modelled with bilinear quadrilateral elements), every single reinforcement bar (modelled with beam elements), and the interface between them (modelled with linear interface elements). Since the beam is symmetric, only half of it was modelled. Furthermore, both support and loading platens were mod- At the symmetry line, the crack band width of the concrete elements was doubled, in order to preserve the fracture energy. This operation rendered another value of the meshdependent parameter in the modified Mazars model, i.e., ε f = 1.17 × 10 −3 . In order to assess the influence of the presence of bond-slip on both the global and local results, a fully-resolved model without interface elements was also constructed. In this model, the reinforcement is assumed to be perfectly bonded with the concrete.
The nonlinear analysis was carried out under displacement control, where the displacement under the loading platen was increased in 100 steps of 0.1 mm. For equilibrium iterations, secant stiffness along with a convergence criterion on the unbalanced forces were used.

Two-scale model
At the large-scale, 6-node triangular elements with 4 Gauss points were used, see "Appendix B" for implementation details. Just as for the single-scale analysis, only half of the beam was modelled, and both the support and loading platens were simulated by tying the nodes appropriately. In view of the fact that the beam was uniformly reinforced, a single type of a RVE was enough to represent the subscale composition. At the symmetry line, the degrees of freedom corresponding to the horizontal displacement and reinforcement slip were locked. Modelling methods and material models described in Sect. 4.1 were employed. A typical two-scale model is schematically shown in Fig. 8.
The size of the unit cell is an important parameter in multiscale modelling, hence 1 × 1, 2 × 2, and 3 × 3 RVEs were used in different FE 2 analyses. Furthermore, to study the sensitivity of the model to macroscopic mesh size, five different meshes were considered, see Fig. 9.
As previously alluded to in Sect. 4.2, the relation between the macroscopic and resolved (subscale) slip is strongly affected by the size of the RVEs. Furthermore, the size of  the macroscopic mesh can also be expected to play a significant role in the analysis. To elaborate on this effect, the macroscopic mesh diameter (average diameter of the area corresponding to a single Gauss point) to unit cell size ratios, h/L , are presented in Table 3.
Similarly to the single-scale analysis, the two-scale simulations were carried out under displacement control, with the displacement under the loading platen increasing in 100 steps of 0.1 mm. In the large-scale analysis, linear stiffness was used for equilibrium iteration. For the subscale analyses, either secant or linear stiffness was used for equilibrium iterations, depending on the RVE. The convergence criterion was set on the unbalanced forces, and if sporadically not met within a predefined number of global iterations, the solution was accepted in spite of lying outside the specified tolerance. On average (in the single steps concerned), the unbalanced force residual was never larger than 1.25 times the intended one. As it was in the same order of magnitude, this was deemed sufficient.

Fig. 11
External load-mid-span deflection relation for two-scale analyses with (ū,s) and without (ū) macroscopic slip variable using 2×2 RVE, compared to fully-resolved (FR) analyses two-scale formulation without macroscopic slip are shown for reference.
As already alluded to in [32], the DD boundary conditions constitute an upper bound on the structural response for the considered problem. This behaviour is noticed also when the macroscopic slip is incorporated. From the results, it can be educed that the incorporation of reinforcement slip as a macroscopic variable has little effect overall on the load- External load-mid-span deflection relation for two-scale analyses with (ū,s) and without (ū) macroscopic slip variable using 3×3 RVE, compared to fully-resolved (FR) analyses deflection response. Even though the results obtained with the smallest unit cell are closer to the fully-resolved solution when the slip is a macroscopic variable, the same cannot be said about the results obtained with 2×2 and 3×3 RVEs. One important cause for the discrepancies between the two-scale and fully-resolved results is the use of Dirichlet boundary conditions on the RVEs. Moreover, it is important to note that even the finest macroscopic mesh is still very coarse compared to the fully-resolved model. Thus, not even the two-scale analysis with the densest mesh can be expected to describe the individual peaks corresponding to localised cracking appearing in the fully-resolved analysis. Furthermore, it can be inferred that the load-deflection results from the two-scale models are not significantly mesh sensitive, since the curves are very near each other for all the meshes. Regarding the effect of the steel-concrete slip in the fully-resolved model, it can be said that the resulting forcedeflection curve is smoother for the case with perfect bond than for the reference case with interface elements. This signals a potentially different fracture development in the model. However, after the cracks have grown, the maximum load is comparable for both fully-resolved models.

Macroscopic strain field
Durability of reinforced concrete structures is heavily influenced by the evolving crack width. Thus, in structural engineering this result is often of main concern, especially in the Serviceability Limit State (SLS) [5], in which the structure can be expected to function most of its lifetime. In SLS, the reinforced concrete structure is usually cracked, but the steel is not yielding. Since in the used constitutive model, the crack opening can be computed by multiplying the strain with the crack band width, it is of interest to study the strain field at SLS. When macroscopic slip was disregarded, the average strain field was smooth and regular, cf. Fig. 13. In the figure, the first principal strain, ε I at load step 62 (resembling SLS) has been plotted for the large-scale meshes. Note that neither large gradients nor strain localisation patterns commonly observed in reinforced concrete structures can be observed at the large-scale. Rather, the strain is "smeared" over a considerably large region.
The incorporation of the macroscopic slip variable had a pronounced effect on the regularity of the macroscopic strain field. The principal strain field,ε I , at load step 62 has been plotted in Fig. 14. In the figure, larger strains and a clear localisation pattern can be seen. In contrast to the two-scale formulation ignoring macroscopic slip, the character of the strain field depended not only on the size of the unit cell, but also on the size of the macroscopic mesh, cf. Sect. 5.3.4 for details. This difference in strain fields was mostly pronounced for finer large-scale meshes.

Crack widths
In the smeared crack approach, the crack widths can be obtained directly by multiplying the strain with the width of the element. In order to study the capability of the developed two-scale model to simulate the evolving crack widths, the comparison with fully-resolved solution is made. To this end, the maximum crack width in SLS was sought in the singlescale solution. Fig. 15 presents the distribution of the first principal strain at load step 62 (resembling SLS),ε I , in the fully-resolved model. The location of the widest crack is indicated by the dot. A corresponding plot for the fully-resolved model with perfect bond is presented in Fig. 16. From the plots, it can be seen that absence of interface elements results in a more distributed cracking pattern with lower principal strains compared to the model with interface elements. To account for distributed cracking in the two-scale models, the integration point with the largest first principal strain, ε I , was sought. For convenience, the search was restricted to the area within a box bounding the process zone in the fully-resolved model, see Fig. 15. Red dots in Figs. 13 and 14 represent the results of the search. At these locations, the loading history was extracted, depending on the type of two-scale model. For the two-scale model ignoring reinforcement slip at macroscale, the strain history was extracted. In the case of the developed two-scale model which allows for macroscopic slip, strain, slip and slip gradient history was extracted. The extracted histories were then imposed on the unit cells via DD boundary conditions. In the ensuing subscale analyses, the maximum crack width within the unit cell 1×1 RVE n el = 23;ū ,s n el = 23;ū n el = 50;ū ,s n el = 50;ū n el = 74;ū ,s n el = 74;ū n el = 108;ū ,s n el = 108;ū n el = 226;ū ,s n el = 226;ū FR FR (perfect bond) Remark The procedure described above might not be necessary, since in theory the results of the subscale analyses for all integration are already available as the output of the FE 2 algorithm. However, it might not be feasible from the 3×3 RVE n el = 23;ū ,s n el = 23;ū n el = 50;ū ,s n el = 50;ū n el = 74;ū ,s n el = 74;ū n el = 108;ū ,s n el = 108;ū n el = 226;ū ,s n el = 226;ū FR FR (perfect bond)

Fig. 19
Maximum crack width, w max , versus mid-span deflection, δ, for the fully-resolved (FR) analyses, and for subscale analyses with 3×3 RVE on different macroscopic meshes. (ū,s) -imposed strain-slip-slip gradient history; (ū) -imposed strain history practical point of view to request output from all integration points (especially for larger unit cells or denser macroscopic meshes). Since the location of the largest principal strain is not known a priori, it is more feasible, from the practical point of view, to find it from the large-scale results in a postprocessing manner; the subscale analysis for a specific unit cell can then be easily repeated.
As shown in [32], the crack widths were underestimated for the DD boundary condition in the original version of the two-scale model. From Fig. 17, it is evident that the crack widths were underestimated for the smallest RVE, for both two-scale models. Moreover, no significant difference could be observed between crack widths computed by the two-scale models. However, a pronounced difference can be observed for 2 × 2 and 3 × 3 unit cells, shown in Figs. 18 and 19, respectively. There, the evolving crack width computed by the two-scale model (with macroscopic slip) agreed much better with the fully-resolved solution, especially for denser macroscopic meshes. It is noteworthy, that for the prescribed strain-slip-slip gradient history, more cracks developed within the RVEs. This is believed to be caused by larger strains computed by the model, as a result of the localised character of the macroscopic strain field (Fig. 14). Furthermore, it can be noticed that the fullyresolved model with perfect bond produced a smaller crack width. This confirms that the steel-concrete bond-slip has a large impact on local results, e.g., crack widths or crack To sum up, incorporation of the macroscopic slip variable had a pronounced effect on the macroscopic strain field. In turn, imposition of the strain, slip and its gradient on the RVE with DD boundary conditions resulted in overall larger crack width, which agreed quite well with the single-scale solution for larger unit cells and denser macroscopic meshes. For the smallest RVEs, the crack width, although slightly increased, was still underestimated.

Macroscopic slip field
Another macroscopic result is the reinforcement slip field, which is plotted for different RVE sizes and large-scale meshes in Fig. 20. It is noteworthy that the regularity of the macroscopic slip field is highly dependent on the size of both the RVE and the large-scale mesh. For coarse macroscopic mesh and small RVEs, the results resembled rather a numerical artefact. For finer large-scale meshes and larger RVEs the field became more regular, i.e., it is possible to see the link with the macroscopic strain field in Fig. 14

Interpretation of scale mixing
In the presented results, we note that the size of the RVE and the resolution of the macroscale mesh have a strong effect on the results. This is a result stemming from incomplete separation of scales; scale-mixing. Since not only the gradient of a macroscopic variable (ḡ), but also the variable itself (s) are imposed on the subscale problem, the results will depend on the size of the underlying RVE. This phenomenon was observed for transient problems in [31], where it was shown that for a given macroscale mesh size, there exists an optimal size of the RVE, which results in the best possible fit to the fully-resolved solution. Also in the present study, there is a strong coupling between the size of the RVE and the macroscale mesh resolution. Note that similar h/L ratio is obtained for n el = 226 (1×1 RVE), n el = 50 (2×2 RVE), and n el = 23 (3×3 RVE) in Fig. 20.
In order to capture the fluctuating fields, a coarse largescale mesh needs a larger RVE, where the fluctuations can develop. On the other hand, the absence of fully developed fluctuation field in a smaller RVE can be compensated by a fine macroscopic mesh, allowing for fluctuations on the higher scale, see Fig. 21  Following the arguments in [31], the resulting accuracy of an approximation using a chosen RVE size and a given macroscopic mesh discretisation, h, can be analysed by studying the two errors 1 stemming from macroscopic FE discretisation and the homogenisation itself. The latter error is the model error introduced by replacing the original problem in Eqs. (11)-(14) by the two-scale counterpart described by Eqs. (42)-(43). When the macroscale fields are sufficiently resolved (i.e., as h → 0), the remaining modelling error comes from the averaging operations, cf. Eq. (30). Clearly this error will increase with the size of the RVE. The discretisation error, on the other hand, is the classical finite element error that is (for a given RVE size) expected to vanish as h → 0.
The purpose of a two-scale formulation is that the modelling error and the discretisation error should cancel out. As illustrated in Fig. 22, an observed quantity Q converges with macroscopic mesh refinement for various sizes of the RVE towards the asymptotic value where only the modelling error prevails. At the limit h → 0, it is thus beneficial to use the smallest possible size of the RVE. However, for coarser meshes (often of engineering interest) the discretisation error may cause a larger RVE to provide a better approximation To summarise, incorporation of macroscopic reinforcement slip results in scale mixing; a combination of the macroscopic mesh size and the unit cell size must be considered in order to obtain the optimal result. The sampling ratio, h/L , was shown to be an important parameter in twoscale analyses, governing regularity of macroscopic fields and accuracy of the results.

Conclusions
In this paper, the framework of Variationally Consistent Homogenisation was used to construct a two-scale model of reinforced concrete, whereby the steel reinforcement, the concrete matrix and the steel/concrete interface are studied in detail. It is then assumed that the reinforcement slip, in addition to the concrete displacements, possesses both a macroscopic and a fluctuation component. The pertinent large-scale problem on the "effective" single-phase solid was outlined in terms of finding the unknown macroscopic displacement and reinforcement slip fields. Within the RVE, the standard linear variation of the macroscopic fields was considered, which corresponds to first-order computational homogenisation. Dirichlet-Dirichlet boundary conditions were considered for the subscale problem on a Representative Volume Element (RVE). A procedure for computing the effective work conjugates, i.e., effective stress, effective transfer stress, and effective reinforcement stress from the subscale problem was outlined.
The subscale response of the RVE subjected to macroscopic reinforcement slip was investigated in analyses of pull-out tests. It was shown that for small RVEs, i.e., with short rebars, the bond-slip law of the interface can be recovered directly from the effective bond-slip relation, since the bond stress is approximately constant along the bar. For larger RVEs, i.e., with longer rebars, the bond stress varies along the bars and the effective bond-slip relation for the RVE does not directly translate into bond-slip law of the steel/concrete interface.
The incorporation of an independent macroscopic slip variable resulted in slightly more involved large-scale problem. The Dirichlet-Dirichlet (DD) boundary conditions constituted an approximate upper bound on the structural response, and the dependency of the result on the size of the macroscopic mesh was negligible. Even though the enrichment with macroscopic slip had little influence on the effective load-deflection relations for a reinforced concrete deep beam, it had a pronounced effect on the character of the macroscopic strain field. When the macroscopic slip is ignored, smooth strain fields "smeared" over the process zone are produced by the two-scale model. With the macroscopic enrichment strain localisation was exhibited by the strain field to some extent. In turn, this localised character resulted in larger strains; these increased the evolving crack width significantly. Although still underestimated by the smallest RVE, the crack width obtained with larger RVEs agreed with the single-scale solution reasonably well.
The macroscopic slip field demonstrated the significance of both the size of the unit cell and the size of the largescale mesh, for the presented problem. A coarse macroscopic mesh with a large RVE, and a fine macroscopic mesh with a small RVE demarcate the potential range, wherein the most meaningful results can be expected. It is noteworthy that the interpretation of macroscopic slip is governed by the size of the RVE.
Apart from verification of the multiscale modelling methodology, the main benefit of multiscale formulations is the potential to save computational effort without compromising the quality of local results. For larger reinforced concrete structures, like e.g., bridges, it may not be feasible to use full resolution in modelling; however, local results such as crack width are still of large importance.
For future work, additional types of boundary conditions could be considered for the RVE problem to ensure more uniform fluctuation of slip and bond stress along longer reinforcement bars. Additionally, three-dimensional reinforced concrete RVE should also be studied in order to be able to model a greater range of large reinforced concrete structures. A first step in this direction can be to extend the the two-scale formulation from 2D solid elements to beam and plate/shell elements.

A.2 Effective transfer stress
The effective transfer stress,τ b , represents the total force transferred from steel to concrete through the interface divided by the area of the RVE, i.e., The transfer force, F Γ is the total force coming from bond stresses distributed along the bars, see Fig. 24. For a twodimensional RVE with one reinforcement bar in x and y direction, the transfer force can be computed as the integral of the bond stresses along the bar, i.e., From the transfer stress, the average bond stress in the bars can be easily recovered as where n i ,S si , andL i are the number, the average circumference and the average length of the reinforcement bars in the i-th direction.

A.3 Effective reinforcement stress
The effective reinforcement stress can be interpreted as the membrane stress for the reinforcement bars only. An analogy between the effective reinforcement stress,σ s , and the effective stress in the reinforced concrete,σ is shown in Fig. 25. The average stress per reinforcement bar (in i-th direction) can be easily extracted from the effective reinforcement stress as σ bar,i =σ s,ii L n iĀsi , (A. 5) where n i , andĀ si are the number and average cross-sectional area of the rebars in i-th direction, and L is the length of the RVE in the direction perpendicular to i.

Appendix B: Large-scale finite element
The fieldsū(x) ands(x) are approximated by shape functions, i.e., ū = N a a,s = N   In the case of plane stress and same order approximation of both displacement and slip fields, the shape function matrix N a = N b = N and the derivative matrices B a , B b can be expressed as

B.1 Internal force vector
Using the finite element approximation in Eqs. (42)-(43), we can compute the internal force vector for the element. Splitting the internal force vector into parts pertinent to displacement and slip fields, we have:

B.2 Stiffness matrix
The stiffness matrix can be computed from linearisation of the finite element equations. Note that slip,s, and its gradient, g, do not depend on the displacement vector a, i.e., ∂s/∂a = ∂ḡ/∂a = 0. Similarly, the strain does not depend on the slip field, i.e., ∂ε/∂b = 0. Linearising the finite element equations, we have