Upscaling of three-dimensional reinforced concrete representative volume elements to effective beam and plate models International Journal of Solids and Structures

Two-scale models for reinforced concrete, where the large-scale problems are deﬁned in terms of Euler– Bernoulli beam and Kirchhoff–Love plate models, are constructed. The subscale problem on the Representative Volume Element (RVE) is correspondingly outlined as ﬁnding the response of the three- dimensional RVE comprising plain concrete continuum, reinforcement bars and the bond between them. The boundary region of the periodic mesh is modelled with special solid elements, which allow for pre- scribing the macroscopic input via strongly periodic boundary conditions in an effective way. The effective response of the reinforced concrete RVEs of different sizes subjected to tension and pure bending is investigated for both effective beam and plate models. A series of experiments on reinforced concrete panels subjected to bending and membrane loads is simulated, and the effective moment–curvature response is studied. Within the developed framework, an arbitrary macroscopic loading in terms of membrane strains and curvatures can be prescribed on the RVE, and the corresponding effective response is obtained, making the proposed formulation feasible for future use in an FE 2 scheme.


Introduction
Concrete is the most widely used construction material in the world, and the use of it is forecast to increase in the future. It is therefore necessary to be able to design and build optimised, safe and sustainable structures. Finite element analysis of concrete structures is a useful tool when developing new structural systems, and it is thus necessary to be able to create robust models capable of providing accurate results efficiently. Reinforced concrete, being ubiquitous in construction around the world typically has cracks appearing already at moderate load levels. Wide cracks can be detrimental to structural performance and thus need to be accurately predicted. This is particularly true for large-scale reinforced concrete structures such as frames, bridges or nuclear reactor containments. In analysis and design of such structures, many simplifications are often made in order to make the computations tractable and to get results within a reasonable amount of time. Therefore, shell and beam elements are often used to approximate the real geometry as they reduce the computational complexity.
Strain localisation in concrete is a complicated phenomenon, and it is not straightforward to obtain accurate predictions. Cracking of plain concrete has been extensively studied in the literature (Hillerborg et al., 1976;Bažant and Oh, 1983;Peerlings et al., 1996;Mazars and Pijaudier-Cabot, 1989;Jirásek, 1998;Grassl and Jirásek, 2006), and there is today a multitude of techniques that treat the cracking problem and provide good quality results (de Borst et al., 2004;Lotfi and Shing, 1995;Ulfkjaer et al., 1995;Grassl and Jirásek, 2010). Computational models, often based on strong (Ottosen and Ristinmaa, 2013;Mosler and Meschke, 2003;Ibrahimbegovic and Melnyk, 2007) or weak (Bažant and Oh, 1983;de Borst, 2002;Rots et al., 1985;Wriggers and Moftah, 2006) discontinuity approaches have made their way to most of the commercial finite analysis codes, and are today widely used by engineers and researchers. These models are usually formulated for two-or three-dimensional solids, but not for beams and shells.
Since concrete structures are almost always reinforced, the situation is further complicated as cracking of the quasi-brittle concrete matrix does not mean structural failure. Instead, structural interaction effects such as tension stiffening (Ben Romdhane and Ulm, 2002;Lundgren and Gylltoft, 2000;Ibrahimbegovic et al., 2010) caused by the bond in the steel-concrete interface must be carefully considered when modelling crack growth. Information about the amount and arrangement of reinforcement is thus vital https://doi.org/10.1016/j.ijsolstr.2020.07.006 0020-7683/Ó 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). and must be accurately represented in the model. All of these reasons make accurate modelling of large-scale reinforced concrete structures challenging and possibly computationally expensive.
In order to tackle both of the described issues, computational multiscale modelling can be employed. More specifically, the FE 2 approach, introduced by Feyel and Chaboche (2000), stands out as an attractive option. In the method, the problem is split into modelling the structure itself at the coarse scale, and the material at the fine scale. Computations are performed on the so-called Representative Volume Elements (RVEs, cf. Gitman et al., 2007), located at the Gauss point level, which are used to model a representative sample of the material at an arbitrary level of detail. The macroscopic loading information is then defined on each RVE via suitable boundary conditions. The response of the RVE is homogenised, whereafter the result is sent back to the coarse scale. It is noteworthy that the FE 2 method addresses both the material model complexity and the potentially large computational cost of analysing large reinforced concrete structures. Within this computational scheme, it is natural to combine the complex fine-scale material behaviour with a large-scale effective model. Furthermore, this approach is very well suited for parallel computing as all RVE problems are independent and, therefore, can be solved simultaneously.
Even though there exists a wide spectrum of multiscale methods and many of them have already been applied to concrete structures (Unger and Eckardt, 2011;Rodrigues et al., 2018;Nilenius et al., 2015;Sun and Li, 2015;Nguyen et al., 2012) or reinforced concrete (Huguet et al., 2017;Sun and Li, 2016;Le et al., 2015;Carbone and Codegone, 2000;Nader et al., 2017) structures, upscaling of the fine scale response of reinforced concrete to effective beam and plate models has not been treated in the literature to the same extent. Studies addressing these problems have focused on different materials, such as masonry (Mercatoris and Massart, 2011;Petracca et al., 2017) or composite structures (Herwig and Wagner, 2018;Coenen et al., 2010;Reinaldo Goncalves et al., 2016;Karttunen et al., 2019). For reinforced concrete, two-scale modelling of frame structures considering effective beam models was studied by Mata et al., 2008. Recently, solid, waffle and hollow core reinforced concrete slabs were analysed in a multiscale setting by Moyeda and Fish, 2019. Still, the FE 2 technique has rarely, to the knowledge of the authors, been used to analyse reinforced concrete structures with large-scale beam and plate models. Studies conducted by the authors (Sciegaj et al., 2018;Sciegaj et al., 2019) considered only twodimensional RVEs, which were upscaled to effective solid elements.
In this work, the first step is taken towards applying the FE 2 procedure to large-scale reinforced concrete structures, which can be modelled with beam and plate elements. To this end, the response of three-dimensional reinforced concrete RVEs is upscaled to effective Euler-Bernoulli and Kirchhoff-Love plate models. First, the developed two-scale model for reinforced concrete (Sciegaj et al., 2018;Sciegaj et al., 2019) is extended within the computational homogenisation framework to cases of having an Euler-Bernoulli beam and Kirchhoff-Love plate models at the macroscale. Next, special boundary elements for prescribing the effective fields on the three-dimensional RVE via periodic boundary conditions are developed. The focus of the paper is placed on the effective response of the three-dimensional reinforced concrete RVE comprising plain concrete continuum, reinforcement bars and the bond between them. Within the developed framework, an arbitrary loading in terms of membrane strains and curvatures can be prescribed on the RVE, and the corresponding effective response is obtained in a straightforward way, making the proposed formulation feasible for future use in an FE 2 scheme.
The remainder of the paper is structured as follows: The singlescale problem is briefly recapped in Section 2. In Section 3, variationally consistent homogenisation is used to derive the twoscale formulation and the necessary scale transitions. In Section 4, the method of imposing the macroscopic fields on the RVE via periodic boundary conditions with help of special boundary elements is described in detail. The modified boundary elements are introduced, and the procedure of generating a periodic mesh for an arbitrary reinforced concrete RVE is outlined. In Section 5, the effective response of the three-dimensional RVEs is investigated in a series of numerical examples comprising both uniaxial tension and bending. Furthermore, the effective moment-curvature response of the RVE is compared to experiments on reinforced concrete panels subjected to bending and membrane loads. The paper is concluded with Section 6, which contains some final remarks and an outlook to future work.

Single-scale problem
The single-scale boundary value problem for a general threedimensional reinforced concrete structure is here briefly recapped, cf. Sciegaj et al. (2018) for a more detailed derivation. The problem domain, X, occupied by the reinforced concrete structure can be split into the concrete part, X c , and the reinforcement part, C int .
The domain boundary can be separated into the essential (C u ) and natural (C t ) parts, with either displacement or traction defined, respectively. A representation of the problem domain along with its components can be seen in Fig. 1.
Denoting the body forceb and the stress in the concrete r c , the momentum equilibrium can be stated as: t :¼ r c Á n ¼t on C t : To retain generality, both loads along and across the reinforcement are considered. Thus, the normal force, N s , is linked to the bond stress, t C , which is distributed along the perimeter of the bar, S s . Even though different reinforcement bars may have different dimensions, it is assumed in the following that the diameter (and consequently the perimeter) are constant for each bar. For the transverse loading in two perpendicular directions, the bending moments, M s ¼ M s;1 M s;2 ½ T , can be expressed as function of the transverse loads k ¼ k ?;1 k ?;2 ½ T . On top of that, it is assumed that there are no end forces applied to the reinforcement bars, i.e. the forces are applied only on the concrete surface. Consequently, the strong form of equilibrium for the reinforcement is: The forces acting on both reinforcement and concrete in the interface are reproduced in Fig. 2. Here, the perpendicular directions represented by the unit vectors e l ; e ?;1 and e ?;2 , are introduced. These vectors can be defined separately for each reinforcement bar comprising the domain C int . For a given reinforcement bar, these directions may be different in different parts of the structure. There is no formal requirement of orthogonality of the reinforcement bars. However, only orthogonal reinforcement grids will be considered in the numerical examples below.
As a next step, we define the displacement fields It is possible to split these fields into longitudinal and transversal components along the interface C int : u s ¼ u s;l þ u s;? ¼ u s;l e l þ u s;?;1 e ?;1 þ u s;?;2 e ?;2 ; ð3Þ u c ¼ u c;l þ u c;? ¼ u c;l e l þ u c;?;1 e ?;1 þ u c;?;2 e ?;2 : For simplicity, we introduce the tensor I ? ¼ I À e l e l ½ , which extracts the transversal part of the displacement vector, i.e., u s;? ¼ I ? Á u s and u c;? ¼ I ? Á u c .
Along the interface, we allow for reinforcement slip, i.e., s :¼ u s;l À u c;l , but the contact deformation in the transverse direction is neglected. Thus, the following interface constraint is adopted: In order to maintain generality, the constitutive relations for concrete, steel, and the interface are considered only implicitly in this work, i.e., they are assumed to be known. The strong forms (1)-(2) and (5) are multiplied by suitable test functions and integrated over the domain. Upon employing the interface equilibrium we obtain the weak forms. Thus, the single-scale quasi-static problem can now be defined as: Find u c ; u s;l ; u s; c dk; u s;? À I ? Á u c ð for suitable trial spaces U c ; U s;l ; U s;? ; L: The following functionals in the system (6)-(9) can then be identified as: (i) Pertinent to concrete: Remark: In the continuous setting, the displacement field in the concrete formally has to be regularised when comparing to the bar/beam elements in Eqs. (7), (8) in order to avoid artificial singularities. One possibility would be to measure the displacement on a rigid cylinder surrounding the reinforcement bars. However, in this study, we restrict to coarse mesh approximations, whereby we assume the finite elements of the concrete to be larger than the physical dimensions of the reinforcement cross-section.

Two-scale problems
In this paper, the term microscale commonly used in multiscale modelling literature has been substituted with the term subscale. This is done to better reflect the chosen level of modelling detail for the RVEs (homogeneous concrete and distinct reinforcement bars can be considered physically macroscopic), and to avoid confusion, as plain concrete has a heterogeneous microstructure, which is not taken into account in this study. Further, the commonly used term macroscale is used interchangeably with the term large-scale in this paper.
The so-called Variationally Consistent Homogenisation (Larsson et al., 2010) is used to derive the two-scale problem. We use the Variational MultiScale (VMS) ansatz and separate the unknown fields u c ; u s;l ; u s;? into the smooth (large-scale) and fluctuating (subscale) parts: with the superscripts M and s denoting the macroscopic and subscale components, respectively. It is assumed that the Lagrange multiplier k lives only at the subscale, and thus the total field is equal to the fluctuation part, k ¼ k s . As the next step, we introduce the running averages, i.e., at each location x 2 X the field is approximated by the volume average on X Ã x ð Þ. Thus, for given functions f X and f C defined on X c and C int , where the running average f Ã is defined as: Utilising the running averages, the fully-resolved problem from Eqs. (6)-(9) is then expressed as Z X a Ã;c ðu c ;du c Þ þ a Ã;l ðu s;l ;u s;? ; du s;l Þ þ a Ã;b ðu s;l ; u s;? ;du s;? Þ þ c Ã ðk; du s; where we introduced the RVE-forms a Ã;l ðu s;l ; u s;? ; du s;l Þ : Consequently, 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. As the integrals are numerically evaluated only at the integration points, RVEs occupy regions X Ã in the vicinity of large-scale integration points, centered around x. Furthermore, it is assumed that the reinforcement bars do not change direction inside the RVE, i.e. the unit vectors e l ¼ e lx e ly e lz Â Ã T and e ?;i ¼ e ?;i;x e ?;i;y e ?;i;z Â Ã T , where i ¼ 1; 2, are constant for each reinforcement bar in the RVE. In order to derive the partial differential equations governing the large-scale problems, as well as the associated scale transitions, we test the fully-resolved problem in (21) with test functions associated with the macroscopic part of the total field, i.e.,

Macroscopic Euler-Bernoulli beam
For the beam model, we introduce model assumptions, presented briefly in Fig. 3. The beam is defined in the xz-plane and spans the length L in x-direction. Its cross-section is assumed to be rectangular with height L z and width L y . The beam is subjected to body loadb ¼bq h i T , where b and q are the longitudinal and transversal loads per length, respectively. At the ends, i.e., at x ¼ 0 and x ¼ L, the beam is subjected to prescribed tractionŝ T . More specifically, we consider a Taylor series expansion of the large-scale fields inside X Ã , centered at x: Upon expanding of (27), we get with L Ã being the length of the RVE in x-direction. Assuming that the applied loads are sufficiently smooth, we discard the higher order terms and express the parts of the right hand side of (35) with help of (31)-(32) as are the prescribed boundary data or reaction forces. The left hand side of (35) is then expanded with help of (31)-(34) to Utilising the symmetry of the stress tensor, we obtain the effective form of (35) as where the effective normal force and bending moment are defined as Finally, the large-scale problem for the macroscopic Euler-Bernoulli beam is defined as: Find with the suitable trial and test spaces.

Remark:
Using the divergence theorem, we have also (assuming the absence of body forces within the RVE) ; where the summation over the discrete forcesR L ;R ? ¼R ?;1R?;2 Fig. 4. These forces together with the boundary tractiont Ã are defined as follows: À M s e ln with the scalar e ln defined as

Macroscopic Kirchhoff-Love plate
For the plate model, we begin with defining the model assumptions, illustrated in Fig. 5. The plate occupies the region A in the xyplane and is assumed to have a uniform cross-section with the height L z . It is subjected to in-plane body loadb p ¼b xby h i T and the out-of-plane loadq. Furthermore, membrane forces N nx ; N ny , shear force V n and bending moments M mx ; M my are acting at a point situated on the external boundary @A with the normal and tangential unit vectors n and m. Similar to body loads, we split the coordinate vector x into the in-plane and out-of-plane components, i.e., The same procedure is applied to the tangential and transversal reinforcement unit vectors, i.e., e l ¼ e lp e lz Â Ã T and e ?;i ¼ e ?;i;p e ?;i;z Â Ã T with the membrane components e lp ¼ e lx e ly Â Ã T and e ?;i;p ¼ e ?;i;x e ?;i;y Â Ã T , respectively. Similarly, we introduce the operator $ p , which denotes the gradient with respect to the in-plane coordinates.
According to the kinematics of this plate model, we have the following relation between the large-scale and subscale displacement fields: In contrast to the beam model, the full VMS ansatz (16)-(18) can be utilised. Following the procedure proposed in (Larsson et al., 2010), we aim at prolonging the large-scale components of the resolved fields from the effective large-scale fields, T . More specifically, we consider a Taylor series expansion of the large-scale fields inside X Ã , centered at Upon expanding (27), we get with A Ã being the area of the RVE projection onto the xy-plane. The right hand side of (53) can be expanded with help of (50) to with the boundary membrane traction, boundary Kirchhoff force, and the boundary moment defined as  Nnx;Nny, shear forceVn, and bending momentsMnx;Mny at a point on the external boundary @A with normal and tangential unit vectors n and m, respectively. N p;n ¼N p Á n;V K n ¼V n þ @ @mM p : n m ½ ;M p;n ¼M p : n n ½ : The boundary membrane forces and bending moments are defined with component representation aŝ Expanding the left hand side of (53) with help of (50)-(52), we obtain: where the effective membrane forces and bending moments are where r p and r z are the in-plane and out-of-plane components of the concrete stress tensor r c , i.e., r p ¼ r xx r xy r yx r yy ! and r z ¼ r xz r yz ! : Finally, the large-scale problem for the macroscopic Kirchhoff-Love plate is defined as: Find with the suitable trial and test spaces.

Remark:
By the divergence theorem, we have also (assuming the absence of body forces within the RVE) wheret Ã;p denotes the in-plane components of the boundary tractionst Ã , and the summation over the discrete forcesR L ;R ? andR M is performed at locations where the reinforcement intersects C Ã , cf.

Subscale problem
In accordance with the VMS ansatz, the fully-resolved problem (21) is tested with functions pertaining to the fluctuation fields. In general, before the boundary conditions on any given RVE are specified, the subscale equilibrium can be expressed in the weak format as with the test spaces defined as follows: Furthermore, in the general case, the boundary terms are given as where the discrete forcesR L ;R ? ;R M , and the tractiont Ã were defined in Fig. 4. Furthermore, the body forces acting within the RVE were neglected. 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. In this paper, we consider periodic boundary conditions, which are imposed on the RVE strongly.
?;1 e ?;i;p " # 4. Periodic RVE problem -finite element modelling In this section, explicit expressions for periodicity constraints are derived for the cases of Euler-Bernoulli beam and Kirchhoff-Love plate models at the macroscale. These constraints serve then as the basis for the finite element implementation, which is also outlined.

Periodicity constraints
In order to introduce periodicity to the subscale unit cell, we divide the RVE boundary into an image part C þ Ã and a mirror part C À Ã . Furthermore, we introduce a mapping u per ðxÞ : C þ Ã ! C À Ã such that all points on C þ Ã can be uniquely associated with corresponding points on C À Ã ; i.e. they have a unique periodic image. It is noteworthy that the definition of the image and mirror parts will vary depending on the macroscale problem. For the beam model, we have a periodic mapping u per;x ðxÞ in one direction. Analogically, a two-directional periodic mapping u per;xy ðxÞ is defined for the case of a plate model at the large-scale. For a visual representation, see Fig. 6.
In terms of scale separation, the periodicity requires that the subscale fluctuations are equal at both image and mirror boundaries, cf. Fig. 7 for a schematic representation for a onedimensional case. Starting from the VMS ansatz (16)-(18), in the first-order computational homogenisation setting we have: ments at the image and mirror boundary of the RVE, respectively. It is noteworthy, that although the equations in this section are stated only for the displacement fields pertaining to concrete, the corresponding formulation holds for the displacement fields in steel reinforcement, and is skipped here for the sake of brevity.

Euler-Bernoulli beam
Keeping in mind the periodicity requirements (u sþ c ¼ u sÀ c and w sþ c ¼ w sÀ c ), we express the total displacement fields at the image and mirror boundary of the RVE and substitute (31)-(34) for the macroscopic parts to obtain the periodicity constraints: where x þ and x À are the x-coordinates of the image and mirror boundary, respectively. It is noteworthy that in the common case of a cuboid RVE, the difference x þ À x À is simply the length of the RVE, L Ã;x . Note that these periodicity conditions depend only on the effective variables (corresponding to macroscopic strains and curvature) and the size of the unit cell.

Kirchhoff-Love plate
Similarly as before, we express the total displacement fields at the image and mirror boundary of the RVE and substitute (50) for the macroscopic parts. It is noteworthy that the macroscopic fields u; v and w do not depend on z. Therefore, the terms pertinent to z in the gradients and coordinates in (50) disappear. Applying the periodicity requirements (u sþ c ¼ u sÀ c ; v sþ c ¼ v sÀ c and w sþ c ¼ w sÀ c ) and dropping the tensor notation, we obtain the periodicity constraints:  where x þ ; y þ ; x À and y À are the x-and y-coordinates of the image and mirror boundary, respectively. It is noteworthy that in the common case of a cuboid RVE, the differences x þ À x À and y þ À y À are simply the lengths of the RVE in x-and y-directions. As before, note that these periodicity conditions depend only on the effective variables (corresponding to macroscopic strains and curvatures) and the size of the unit cell.
In the remainder of the section we simplify the notation used for the macroscopic deformation gradient and curvatures tensors as follows: Similarly, we denote the distances between the image and mirror boundaries y þ À y À ¼ L Ã;y :

Periodic mesh generation
The easiest way of satisfying the periodicity constraints (73)-(74) or (75)-(77) in the RVE model, is to apply them strongly to specific finite element nodes. As a result, the finite element mesh has to fulfil certain requirements. First, all nodes on the image boundary must have a periodic counterpart on the mirror boundary. Secondly, for any node on the image boundary, it must be possible to uniquely identify its periodic counterpart. The method presented here is inspired by the works of Yip et al. (2005) and Grassl and Antonelli (2019) and consists of four stages.
First, points are randomly generated within a specified geometry, usually a cuboid. While there are no requirements on the points generated within the region, it is required that as soon as a point is placed on a surface, another point is automatically placed in the same location on the opposite surface. In this way, it is assured that regardless of which boundaries are later considered as image or mirror boundaries, a periodic image exists for any point located on the surface of the unit cell. Secondly, given the list of point coordinates, Delaunay tessellation is performed. In this stage, the solid finite elements, which will later constitute the concrete phase, are generated. In the next stage, the dual Voronoi tessellation is performed for the same point input, and the Voronoi polyhedra associated with each point are generated. Those polyhedra bound regions in space within which it is closer to the associated point than to any other point in the neighbourhood. Although not formally needed to generate the finite element mesh for concrete, this tessellation is later used to include reinforcement and the steel-concrete interface in the RVE, which is done in the last stage. Since the reinforcement bars are modelled with beam elements, the discretisation is rather straight-forward due to the simple geometry of the rebars. Given the start and end point of any reinforcement bar and the Voronoi polyhedra data, the intersection of the rebar line with the polyhedra can be found. Having found the entry and exit point within a specific polyhedron associated with a concrete node , the reinforcement node is placed midway between the intersection points. To model the interface, a node-to-node interface element is placed between the concrete and reinforcement node (see Fig. 8). Lastly, for all elements with nodes on the image boundary, the periodic images of these nodes are found and the information is stored for future use.
It is noteworthy, that this design allows to quickly generate periodic mesh for a reinforced concrete RVE with arbitrary directions and reinforcement layout. The coordinates of the region occupied by the RVE, and the start and end points of each reinforcement bar, is all information needed. Although it is assumed that the reinforcement bars will go throughout the RVE, it is also possible to include rebars which are partly or completely contained, i.e. which start or end, inside the RVE.

Modified boundary elements
In order to resolve the periodicity constraints within a finite element model, special modified boundary elements were implemented in the open source C++ finite element code OOFEM (Patzák, 2012). Formulation of these finite elements, located at the image boundary of the RVE, takes advantage of the fact that the degrees of freedom on the image boundary are not free. Instead, due to the periodicity constraints (73)-(74) and (75)-(77), they depend on the periodic degrees of freedom located at the mirror boundary, the effective variables coming from the macroscale, and the size of the unit cell. Therefore, these degrees of freedom can be eliminated from the equation system during solution, and later easily reconstructed for the purpose of postprocessing. The strategy is represented schematically for a tetrahedral element in Fig. 9, which shows an example of the modified element with three periodic nodes. In the picture, nodes i 0 ; j 0 and k 0 are situated at the image boundary, while the remaining node l 0 is located inside the RVE.
It is noteworthy that in the general case, any number of element nodes can be situated at the image boundary. Also note that the information whether the specific node has a periodic image is readily available from the mesh generation stage. We gather the degrees of freedom of the boundary element in the element vector and the ''master" degrees of freedom in element vector where C is a vector containing the effective variables coming from the large-scale. Contents of this vector depend on the choice of the macroscopic model. For the Euler-Bernoulli beam model we have while the Kirchhoff-Love plate model yields Remark: In addition to the macroscopic variables, rigid body motions must be prevented. This means constraining 3 translations and 2 rotations for the Euler-Bernoulli beam model and 3 translations for the Kirchhoff-Love plate model.
We are now in position to establish the relationship between u 0 e and u e with the help of periodicity constraints (73)-(74) and (75)-(77) as with T e representing the periodicity-related linear transformation matrix. More explicitly, we have where z a is the z-coordinate of node a 0 , and p ax is an indicator variable equal to 1 if node a 0 has a periodic image in x-direction and 0 otherwise. Analogically, p ay is an indicator variable equal to 1 if node a 0 has a periodic image in y-direction and 0 otherwise.
During assembly of the stiffness matrix and internal force vector for the system, the contributions from the boundary elements are computed as where K 0 e and f 0 int are the element stiffness matrix and internal force vectors with respect to the nodes located at the image boundary, and T e is the element transformation matrix. Effectively, all the degrees of freedom at the image boundary are eliminated from the equation system, but they are easily recovered with (82).
Regarding practical implementation aspects, few things are worth mentioning. First of all, an analogous treatment can be applied to beam elements which are used to model reinforcement bar. In this case, the equations are simplified, since we know that at most only one node will be located at the image boundary. The analogous derivation for the beam element is not repeated here for the sake of brevity. Secondly, the macroscopic effective variables are stored as degrees of freedom of a global control node, which is shared among all boundary elements. Therefore, each boundary tetrahedral element is defined by 5 nodes (3 nodes for boundary beam element). Note that the control node is global for a specific RVE problem and not global for the macroscopic problem, i.e. every unit cell has a unique control node which stores the macroscopic data associated with it. Third, an arbitrary effective strain (curvature) history can be imposed on the RVE in an elegant and effective way by prescribing the effective variables as boundary conditions in the control node. Lastly, the effective work conjugates (which translate to effective membrane forces and bending moments) are obtained automatically as reactions in the control node. It is noteworthy that the effective work conjugates can also be imposed on the RVE by prescribing the associated forces, which can be useful in stress control setting at the macroscale.

Numerical examples
In order to ascertain the feasibility of the developed threedimensional reinforced concrete RVE model for use in a potential FE 2 scheme, it is necessary to verify its response under some common loading cases. In this section, the response of RVEs in terms of the effective membrane forces and bending moments as well as subscale fracture patterns is investigated. First, RVEs corresponding to both effective beam and plate models are subjected to uniaxial tension. Next, the RVEs from both effective models are subjected to uniaxial bending. Lastly, a series of experiments on reinforced concrete panels subjected to bending and membrane loads is modelled with the proposed technique, and the effective moment-curvature response is compared. All the simulations presented here were carried out in the open source C++ finite element code OOFEM (Patzák, 2012).

Uniaxial tension and bending for effective beam and plate models
Before introducing the individual tests, the assumptions behind the constitutive models for all materials are briefly presented, as these choices were common across all the simulations in this subsection. For the concrete, an isotropic damage model with and exponential softening in tension was used. The Rankine strength envelope was used, i.e. tensile stresses were limited while a linear elastic response in compression was considered. For the cracking, the smeared crack formulation, with a crack band width equal to element size, was used. The concrete grade was assumed to be C30 and the necessary material parameters were obtained from the Model Code (Fib, 2013): elasticity modulus of 33.6 GPa, Poisson's ratio of 0.2, tensile strength of 2.9 MPa and mode I fracture energy of 140.5 N/m. For the reinforcing steel, the von Mises ideally plastic material model was used with Young's modulus of 200 GPa, Poisson's ratio of 0.3 and the yield strength of 500 MPa. The bond-slip behaviour of steel-concrete interface was modelled according to the recommendation of the Model Code (Fib, 2013), cf. Fig. 10. The interface characteristics s 1 ¼ 1 mm; s 2 ¼ 2 mm; s 3 ¼ 6:5 mm; s bmax ¼ 15:41 MPa and s bf ¼ 6:16 MPa were adopted.

Uniaxial tension
In the first series of tests, uniaxial tension was considered as the macroscopic load case. First, the effective Euler-Bernoulli beam model was considered. An example of the unit cell for this model can be seen in Fig. 11.   The beam had a square cross-section with L y ¼ 0:1 m and L z ¼ 0:1 m and was reinforced with one steel bar (diameter / ¼ 16 mm) in the middle of the cross-section, resulting in the effective depth d ¼ 50 mm. In order to investigate the effect of RVE length on the effective response (and on the crack spacing), eight different unit cells sizes L Ã ¼ L Ã;x were considered. Namely, the RVEs with lengths 0:2 m ðn el ¼ 8034Þ; 0:3 m ðn el ¼ 12154Þ; 0:4 m ðn el ¼ 15956Þ; 0:5 m ðn el ¼ 19567Þ; 0:6 m ðn el ¼ 23657Þ; 0:7 m ðn el ¼ 27912Þ; 0:8 m ðn el ¼ 31595Þ and 0:9 m ðn el ¼ 35368Þ were studied. In the parentheses, the number of finite elements for each unit cell is given.
Uniaxial tension was imposed by prescribing the degree of freedom corresponding to the effective strain E xx in the control node of the RVE model. The simulations were run in displacement control, where the total macroscopic strain of E À xx ¼ 2:75 Â 10 À3 was prescribed in 275 steps. In order to ensure that deformation is purely tensile, the degrees of freedom associated with macroscopic rotation E zx and curvature K xx were set to 0 for both beam and plate models. The effective normal force was then obtained as the corresponding reaction in the control node. The results for the effective Euler-Bernoulli beam model can be seen in  In the graph, the effective normal force was normalised with respect to the force causing yielding of the reinforcement bar, F y ¼ A s f y , where A s is the cross-section area of the reinforcement bar and f y denotes the yield strength of steel. The macroscopic strain was normalised with respect to the yield strain of the reinforcement, e y ¼ f y =E s , where E s is the Young's modulus for steel. Furthermore, the average crack spacing, s m , was estimated for each unit cell and is indicated in the figure. It can be concluded that the Fig. 14. Normalised effective normal force versus normalised macroscopic axial strain for beam RVEs with lengths 0.8 m and 0.9 m. Strain localisation patterns and the average crack spacings (sm) indicated. effective response of the RVEs agrees with the well-known response of reinforced concrete in tension. First, the response is linear elastic until first cracking occurs. After that, the concrete experiences softening and load is carried by the uncracked concrete and reinforcement, which explains the change in slope of the response. If the geometry allows for it, i.e. if the RVE is large enough, more than one crack can form in the model, which is reflected in the strain localisation patterns for each model. After cracks have formed, the load can increase until the stress in the steel reaches the yield stress. Subsequently, the deformations will increase at constant load due to the assumption of ideal plasticity. Qualitatively, it can be concluded that the effective response of the RVEs follows this theoretical development and the strain localisation patterns correspond to a large extent to what can be expected in reality. The exact shape of the force versus strain response depends on the number of cracks, but overall the response is independent of the unit cell after the initial cracking is resolved.
As the crack spacing in a reinforced concrete tie in tension is governed by the geometry, reinforcement layout and material properties, ideally one constant value of crack spacing would be expected for all RVEs. However, in an RVE, only an integer number of cracks can be observed. Therefore, if the length is not large enough for another crack to fit inside, no more cracks will appear in the RVE. In all simulations, the first crack appeared at the boundary, while the first intermediate crack appeared in the middle of the RVE (first in the 0:4 m long RVE), thus dividing the unit cell into two equal parts. Subsequent cracks were observed in the RVE only when the length of those parts was large enough for another intermediate crack to fit inside on both sides, i.e., 3 intermediate cracks were only observed in the 0:8 m and 0:9 m long RVEs. For the shorter RVEs (0:5 m; 0:6 m and 0:7 m long) only one intermediate crack was observed in the middle. This can be seen in Figs. 12-14. From the figures it seems that the total number of cracks is always doubled if the RVE is long enough (either 1, 2 or 4 cracks were observed). As the crack spacing is not known a priori, it might be necessary to use a large (long) RVE to confirm the correct crack spacing.
As a next step, the influence of the element size within the RVE of a fixed length was investigated. To this end, the 0.4m long RVE was modelled with various size of finite elements. The tensile test described earlier was then repeated for all the models. The results for five different finite element meshes can be seen in Fig. 15, where the normalised force versus strain response was plotted for the different element sizes. The finite element models contained 2056; 15956; 31374; 47200, and 124818 elements in total, respectively. The approximate average element sizes were then computed for each model by relating the volume of the RVE to the number of concrete finite elements, assuming that all concrete elements are regular tetrahedra. From the graph, it can be concluded that there is good agreement between the different mesh sizes when it comes to the effective response. Even though the secondary crack appears at different stages of the loading for different mesh sizes, the slopes both at the initial stage and after softening agree with each other. Furthermore, the strains tend to localise in the same way regardless of the mesh size, and thus the average crack spacing remains the same across the meshes.
To summarise, the initial crack formation and the ultimate plateau are independent of the RVE size. While qualitatively the same, the transition between these two states depends on the RVE size (number of fractures) and the mesh size.

Uniaxial bending
In the subsequent test series, uniaxial bending was considered as the macroscopic load case. Previously studied RVEs were used as a starting point for these simulations. The main modification consisted of lowering the position of the reinforcement bar, as this reflected better what is done in practice for reinforced concrete beams. Therefore, the effective depth d was changed from 50 mm  Subsequently, similar RVEs corresponding to large-scale Kirchhoff-Love plate model subjected to uniaxial bending were considered, cf. Fig. 16. The plate had a thickness of L z ¼ 0:1 m and was reinforced with a reinforcement grid consisting of 16mm bars in both directions spaced 0.1 m apart to agree with the previously described beam geometry. The reinforcement grid was placed closer to the bottom of the cross-section, with the effective depths d x ¼ 72mm and d y ¼ 56 mm. In this case, the unit cells were assumed to be cuboids with a square base, i.e., it was assumed that Similarly, three different RVEs with L Ã ¼ 0:2 m; L Ã ¼ 0:4 m and L Ã ¼ 0:8 m were considered for the simulations. The unit cells comprised 16044; 64390 and 255248 elements, respectively.
Uniaxial bending was imposed by prescribing the degree of freedom corresponding to the effective curvature K xx in the control node of the RVE model. The simulations were run in displacement control, where the total macroscopic curvature of K xx ¼ 0:1 m À1 was prescribed in 400 steps. In order to eliminate overconstraining, the degrees of freedom corresponding to macroscopic strains remained free, i.e., the RVEs could extend and contract freely. The rotations E zx for the beam model as well as E zx and E zy for the plate model were set to 0 in order to simulate uniaxial bending. The effective bending moment was then obtained as the corresponding reaction in the control node.
The results for the effective Euler-Bernoulli beam model can be seen in Fig. 17 while the results for the effective Kirchhoff-Love plate model are presented in Fig. 18. In the figures, the effective bending moment was normalised with respect to the ultimate moment, M y , of the cross-section. The moment was calculated approximately as M y ¼ 0:9A s f yd d, where A s and f yd are the reinforcement area and yield strength, respectively. The corresponding curvature at yielding, j y , was used to normalise the macroscopic curvature. The curvature was found in a post-processing step and corresponds to the macroscopic value at the stage when yielding in the reinforcement happens for the first time. Furthermore, average crack spacings, s m , were estimated for the different unit cells and are indicated in the figures. By looking at the graphs, it can be concluded that the response of the RVEs for both macroscopic models agrees with the theoretical response of reinforced concrete section under bending. First, the response is linear elastic until first cracking occurs (stage I). In contrast to axial loading, the cracking starts at the bottom of the cross-section, leaving the concrete at the top uncracked and subjected to compressive stresses. The tensile stresses are after cracking carried by the reinforcement bars and the concrete between the cracks (tension stiffening). Thus, the slope of the effective response changes to reflect the weaker cracked parts in the model (stage II). As the load increases, new cracks form at the bottom of the cross-section and the already existing cracks grow. The load can increase until the reinforcing steel reaches yielding (stage III). At even later stage the section fails either by rupture of the steel, or by crushing of the concrete in the top of the beam. However, this later stage is not described in the model, because of the assumptions of ideal plasticity of the steel and the Rankine strength envelope for the concrete. As a result, the curvature will increase at constant load due to the assumption of perfect plasticity.
Qualitatively, it can be concluded that the effective response of the RVEs under uniaxial bending for both large-scale beam and plate models follows the expected theoretical three-stage development and the strain localisation patterns correspond to what can be expected in reality. There does not seem to be a large dependence on the size of the unit cell after the initial cracking is resolved. For the plate models, the number of cracks formed per length of the RVE was slightly higher for the larger unit cells than in the corresponding beam RVEs. In case of the beam RVEs, the cracks were short and as result could be considered to be straight. In case of the plate RVE, the cracks were no longer straight, and the inter-crack distances varied along the width of the RVE and it was also possible to observe branching of the cracks. However, the average crack spacing showed a rather small variance between the unit cells of different lengths.

Remark:
It is noteworthy that the choice of periodic boundary conditions has an impact on the cracking that can occur within the RVE. For the case of beam, periodicity is enforced only in one direction (at the two opposite faces normal to the longitudinal direction). As such, this does not prevent the formation of inclined cracks which cross the other surfaces of the RVE. The situation is slightly Fig. 17. Normalised effective bending moment versus normalised macroscopic curvature for three different beam RVEs. Strain localisation patterns and the average crack spacings (sm) indicated. different for the case of plates, where the periodicity is twodimensional, and the proposed implementation will produce crack patterns aligned with the (in-plane) coordinate axes. As a result, inclined cracks are generally not possible without special treatment. To alleviate this issue, it is possible to align the periodicity with the direction of the crack, which has been studied e.g., in Coenen et al. (2012), Svenning et al. (2017) and Svenning et al. (2019). Such a treatment might prove necessary if inclined crack growth, e.g., close to the supports of the structure, is to be studied. It is noteworthy that the two-scale method proposed by Svenning et al. (2019) does not require a priori knowledge of the crack direction, as this issue is handled by the solution procedure. Moreover, different points in the large-scale domain can have different periodicity orientations.

Bending and membrane loads -experimental comparison
Although the previously presented virtual simulations agree qualitatively and quantitatively with what can be expected for reinforced concrete, it is useful to study the feasibility of the devel-  oped model in simulating real reinforced concrete sections. A series of experiments suitable for this purpose was conducted at University of Toronto and reported in Polak and Vecchio (1994). In the tests, reinforced concrete panels were subjected to bending and membrane loads. A system of in-plane and out-of-plane actuators applied the loads along the edges of the panel and made it possible to subject the panel to a wide range of combinations of bending moments, in-plane and out-of-plane forces. The layout of the panel tests together with the geometry and reinforcement layout is schematically presented in Fig. 19, based on the report from the study.
Three experimental loading cases, named SM1 through SM3, were analysed, cf. Table 1. In the SM1 experiment, the panel was tested in uniaxial bending. The SM2 test comprised uniaxial bending together with in-plane tension in x-direction as well as inplane compression in y-direction with the loads fulfilling the condition Lastly, the response under biaxial bending was studied in the SM3 test, with the bending moment in x-direction being 3.2 times larger than the corresponding moment in y-direction, i.e. M xx ¼ 3:2 M yy . The tested panels had a square shape and were uniformly reinforced, with identical reinforcement placed on the top and bottom of the cross-section. The thicker rebars oriented along the x-direction were placed closer to the surface. The reinforcement grid had a uniform spacing of 76 mm in both directions. Three-dimensional unit cell with a square base of length 836 mm, marked in Fig. 19, was considered for the simulations. The unit cells contained 11 reinforcement bars per layer in each direction and comprised 110892 elements.
As previously, an isotropic damage model with Rankine strength envelope and exponential softening in tension was used for the concrete. Material parameters for both the concrete and the steel-concrete interface are summarised in Table 1. In the original report, only the compressive peak strength and strain, as well as the tensile strength obtained from standard cylinder tests, were given. Due to lack of information on the additional parameters needed for the constitutive model, they were estimated according to the guidelines in the Model Code (Fib, 2013).
The constitutive behaviour of reinforcing steel was reported as trilinear with perfect plasticity until the hardening strain is reached. After that, a linear hardening regime is entered until the ultimate strain is reached. In the simulations, this behaviour was simplified to bilinear elastoplastic model, with linear hardening starting directly after yielding. The material parameters for the reinforcing steel along with some interface characteristics are summarised in Table 2.
In the simulations, arc-length control was used with the reference loading prescribed in the corresponding degrees of freedom in the control node. In this way, the solver finds a deformation of the model such that the resulting reactions satisfy the given constraints, e.g. moment or force ratios. In the SM1 analysis, a reference moment of 1 Nm was applied to the degree of freedom corresponding to K xx in the control node. In order to introduce the in-plane loads in the SM2 test, reference loads of 4 N and À4 N were applied to the degrees of freedom corresponding to E xx and E yy . In conjunction, a bending moment of 1 Nm applied to the degree of freedom corresponding to K xx . Similarly, in the SM3 simulation, a reference moment of 3.2 Nm was prescribed in the degree of freedom corresponding to K xx and a reference moment of 1 Nm was prescribed in the degree of freedom corresponding to K yy . The load factor, output directly by the solver, can then be interpreted as the macroscopic load. The actual values of the effective curvatures and in-plane strains can be obtained directly from the ''displacements" of the control node.
The results for the SM1 through SM3 simulations are presented in Fig. 20a-c, with the experimental results marked with crosses. Even though the results overestimate the effective bending moments slightly, the stiffnesses (slopes) agree well with experimental results. The slight discrepancy might be caused by the lack of necessary information about the concrete. Furthermore, it is noteworthy that the proposed formulation contains ''extra" concrete in the volume of the unit cell, as the reinforcement is simply added to the model. For the tested cases, this adds around 4% of superficial concrete volume. In view of these facts, it is considered that the RVE model reflected the experimental behaviour rather well.
Even though crack patterns for plate RVEs subjected to uniaxial bending have already been studied in previous sections, it is still insightful to look at the crack pattern and spacing for the case of  biaxial bending, which is simulated in the SM3 test. Fig. 21 shows the crack pattern along the width on the bottom surface of the unit cell at the initial yielding of reinforcement, i.e. at the macroscopic curvature K À xx ¼ 12 Â 10 À3 m À1 . The average crack spacings in xand y-directions are indicated as well. From the figure, it is noteworthy that cracking in both directions is captured. The average spacing between cracks parallel to y-direction was reported to be 100 mm, while the average spacing between cracks parallel to xdirection was reported to be 200 mm. The corresponding values of 167.2 mm and 209 mm were obtained in the simulation, and are considered to reflect the experiments moderately well.
It is noteworthy that the only difference between the simulations (apart from the minor discrepancy between material parameters) was the load vector of the control node. As a result, it can be concluded that an arbitrary loading state (in terms of in-plane loads and bending moments) can be prescribed on the RVE and the corresponding effective response is obtained in a reliable and straightforward way. This can in turn function as an ''ad hoc" constitutive model of reinforced concrete, making the proposed formulation FE 2 friendly.

Conclusions
In this paper, the response of three-dimensional reinforced concrete RVEs was homogenised to macroscopic beam and plate models. To this end, special modified finite elements were constructed. The new RVE formulation makes it possible to apply an arbitrary combination of macroscopic membrane strains and curvatures Fig. 20. Moment-curvature relations for SM1 (a), SM2 (b) and SM3 (c) test simulations. Experimental results from Polak and Vecchio, 1994. via periodic boundary conditions. The performance of the RVEs coming from both effective beam and plate models was investigated in a series of analyses, where the unit cells were subjected to uniaxial tension and bending. It was found that the effective response of the RVE followed the expected response of a reinforced concrete cross-section in tension and bending. Furthermore, the influence of the size of the RVE on the effective response was negligible. Similarly, the effective response and crack patterns were consistently captured by the model for various finite element mesh sizes. Subsequently, the RVE was used to simulate experiments, in which reinforced concrete panels were subjected to bending and membrane loads. The simulations predicted the moment-curvature relations well, with the response being slightly overestimated. Cracking in both directions was captured correctly by the model, and the average crack spacings were reflected in a satisfactory way. Within the developed framework, an arbitrary loading in terms of membrane strains and curvatures can be prescribed on the RVE and the corresponding effective response is obtained in a straightforward way, making the proposed formulation feasible for future use in an FE 2 scheme.
The main benefit of the proposed method is the reduction in the number of degrees of freedom from full three-dimensional reinforced concrete representation considering embedded reinforcement and bond-slip to effective beam and shell models. For the proposed method to be applicable to real-world applications, the size of the macroscopic elements should be much larger than the crack spacing observed in the structure. The intended application of this method are very large reinforced concrete structures which are often modelled with shell elements, such as bridges or nuclear reactor containment buildings. Nevertheless, as the concept is still valid for other materials, the proposed method could also be applied to thin textile reinforced shells, which exhibit much finer cracking networks. Furthermore, the FE 2 models are a necessary basis for reduced order modelling which can improve the numerical efficiency of engineering simulations. The homogenisation techniques can also be used for identification of effective macroscopic beam/shell models. Last, the proposed RVE modelling approach could also be used for investigating parts of the structure where the impact of reinforcement detailing is of interest.
As future work, this subscale RVE model will be integrated in the full FE 2 scheme, making it possible to simulate the behaviour of large structures and at the same time enabling the study of crack growth in detail in the RVEs. Moreover, before the FE 2 scheme is applied to analyses of real structures, it is necessary to investigate the validity of the scale separation assumption for this particular application. In the case of reinforced concrete, scale separation is rather small. It might therefore prove difficult to model some more complex structures, where e.g. high strain gradients are present already at the macroscopic level. In this setting, it is also important to extend the current formulation to effective Timoshenko beam and Mindlin plate models to allow for modelling a wider range of structures, for which the assumptions of Euler-Bernoulli or Kirchhoff-Love models are too restrictive. Furthermore, in the current model, the reinforcement slip varied only locally, i.e., at the RVE level. It was shown in previous work by the authors (Sciegaj et al., 2019) that it is necessary to allow for reinforcement slip variation at the macroscale. A similar strategy can be applied also to the effective beam and plate models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.