On the Calibration of a Numerical Model for Concrete-to-Concrete Interface

The study was devoted to the numerical modelling of concrete-to-concrete interfaces. Such an interface can be found in many modern composite structures, so proper characterisation of its behaviour is of great importance. A strategy for calibration of a model based on cohesive finite elements and the elastic-damage traction–separation constitutive law available by default in the Abaqus code was proposed. Moreover, the default interface material model was enhanced with the user-field-variables subroutine to include a real strength envelope for such interfaces. Afterwards, the modelling approach was validated with numerical simulation of the most popular tests for determining the strength characteristics of concrete-to-concrete interfaces: three-point bending beam with a notch, splitting bi-material cubic specimens, and slant-shear tests. The results of own pilot studies were used as well as those reported by other researchers. The performed simulations proved the accuracy of the proposed modelling strategy (the mean ratio of ultimate forces obtained with numerical models and from experiments was equal to 1.01). Furthermore, the presented examples allowed us to better understand the basic test methods for concrete interfaces and the observed mechanisms of failure during them.


Introduction
Numerical modelling of concrete, cf. [1][2][3], and its contact zones with other buildings material is and will continue to be a current topic due to advances in concrete technology and material engineering [4]. In every reinforced concrete member, there are steel-toconcrete interfaces, which were thoroughly studied in reference [5]. The interfaces could also be found in composite structures as tubular columns filled with concrete [6] or in structures strengthened with composite sheets [7]. The last mentioned reference shows how effective a modern finite element (FE) system can be for structure assessment, since the failure modes and force-deflection curves for complex retrofitted beams were perfectly captured in this study.
The present study puts forward concrete-to concrete interfaces. Interfaces of different types of concretes, whose strength properties vary significantly, can often be found in modern structures, cf. Figure 1. For example, they occur in composite structures as slimfloors [8,9] or composite T-beams [10]. In the former, they influence the structure load capacity and deformability tremendously, since they are located at the interface between prefabricated hollow slabs, hybrid beams, and cast in situ overlay concrete, so their total area is of great value. Such types of interfaces can also be found in high buildings at the connections between columns and floors [11]. The columns are usually made of highperformance concrete (HPC), whereas the concrete class of floor can be remarkably lower. A concrete-to-concrete interface is also formed as a result of structure strengthening by casting additional concrete layers [12]. Concrete-to-concrete interfaces usually have lower strength properties than the concrete layers that form it. Therefore, the proper and robust numerical model for the interface is crucial to simulate the behaviour of many composite structures. It has to precisely predict the strength of the interface (in a complex traction stress state, i.e., with normal and shear tractions) as well as its post-cracking behaviour. Such models are implemented in popular FE codes like ATENA [13], ANSYS [14], or Abaqus [15]. The calibration of an interface model in the first code, which is dedicated to the modelling of concrete structures, is easier than in the other two general-purpose FE systems. Especially, the standard damage initiation criterion in the ANSYS or Abaqus code is not suitable for concrete under a complex traction state.
To determine the mechanical parameters of the interface model, the results of different strength tests can be used. They were widely described in many references [12,[16][17][18][19].
In general, such tests aim to induce a homogenous traction stress state at the interface: tension or shear dominated (with compressive normal traction stresses as well). Schemes of the most popular test methods are shown in Figure 2. A pull-off test is often used in a quick assessment of bond strength [20] and can be used for in situ conditions as well. Modifications of standard laboratory tests for homogenous concrete like splitting or direct tension are also used for samples prepared in a laboratory [18,21]. In Figure 2, the most popular tests for determination of interface shear strength were presented [22,23], but new testing methods have also been developed [24]. The easiest way of testing the interface strength under compression and shear is a slant-shear test [25]. To determine the fracture energy of an interface, methods initially developed for homogenous concrete are used, like the three-point bending sample with a notch [26,27] or the wedge splitting test [28]. For a quick assessment of the strength properties of an interface, it is convenient to introduce the bond efficiency coefficient α int [16,29] which can be calculated as follows: where: f int and f w are the strength of the interface and the weaker concrete layer obtained from the same test method, respectively. The bond strength between layers of concrete cast at different times is affected by many factors described in numerous references, cf. [16,25,30,31]. After the literature survey, the following factors can be considered as the most important: properties of the overlay and substrate concretes, surface characteristics (its cleanliness and mechanical preparation), compaction of concrete mix, and curing conditions, which influence shrinkage strains. The influence of the shrinkage is hard to capture with specimens of small dimensions, which was clearly presented in reference [21]. There are a few recent studies concerning the calibration of a numerical model for the concrete-to-concrete interface. Frenzel and Cubarch [32] performed shear tests of a concrete interface with regular, infra-lightweight, and foam concrete. On this basis, they calibrated a numerical model for a concrete-to-concrete interface subjected to pure shear in the ATENA code and obtained satisfactory agreement between the simulation and test results. Farzad et al. [33] conducted three types of tests: flexural, direct shear, and slantshear ones, which enabled them to obtain strength characteristics of ordinary-to-ultra-high performance concrete (UHPC) interface in different traction stress states. Afterwards, they proposed parameters for the concrete damage plasticity model (CDP) in the Abaqus code for fictitious continuum material corresponding to the interface behaviour. However, such an approach is not recommended by the Abaqus manual [15]. Using continuum material models with an interface element is appropriate in the case of flexible interfaces (with lower stiffness than the stiffness of connected parts) and when it is possible to determine the strength and deformability of the material, which constitutes an interface in laboratory tests. Such a situation does not occur for a concrete-to-concrete interface-this kind of connection is obviously initially rigid, and it is impossible to extract its material due to its infinitesimal thickness. According to the manual [15], a material model based on a traction-separation relation should be used for such kinds of interfaces. Valikhani et al. [22] calibrated a model based on cohesive elements and the traction-separation law available in ATENA for a concrete-to-UHPC interface. Only shear loading was considered in this study. The latest study by Yuan et al. [34] concerns a rock-to-concrete interface. It contains detailed research on composite beams made of concrete and rock as well as the proposition of calibration of the traction-separation law for interface subjected to tension (mode I of fracture) implemented in the ANSYS software [14]. All in all, the authors of this study could not find in the literature a complete strategy for calibration of a concrete-to-concrete interface subjected to various loading conditions (combination of normal and tangential traction stresses), especially using cohesive elements and the traction-separation relation available by default in the Abaqus or ANSYS code.
The main aim of this study was to propose a complete strategy for calibration of a concrete-to-concrete interface according to the mentioned assumptions. Simulations were performed with the Abaqus code [15] using cohesive elements and the traction-separationtype material model. The standard constitutive model was enhanced with a simple USDFLD user subroutine coded in FORTRAN, which facilitates the correct prediction of the interface strength subjected to compressive normal forces. The issue of fracture energy in two modes of fracture (tensile and shear) was discussed and implemented in the analyses. The results of experimental tests found in the literature (i.e., in references [21,25,26]), as well as the results obtained during own pilot studies, were used to validate the proposed strategy. The validation part was divided into four case studies (CS). The following types of concreteto-concrete interface laboratory tests were examined: three-point bending of a beam with a notch, splitting of cubic specimens, and slant-shear tests on prisms.
The major novel elements presented within the article are as follows: a complete strategy for calibration of interface traction-separation laws for concrete-to-concrete interfaces available by default in the Abaqus code and the detailed non-linear simulation of the most popular strength tests of concrete-to-concrete interfaces. They took into account: advanced constitutive models for concretes and their interfaces as well as detailed geometry of the loading devices. It is worth mentioning that the progress of delamination of bi-material specimens is hard to follow with modern optical measurement systems due to their brittle failure modes, so only numerical studies could give a wider insight into this process at this point [23,27].

Notched Beam by Chen et al. (2021)
Chen et al. [26] performed a three-point bending test on high-strength and lowstrength concrete and two types of their interface, namely, without special surface preparation and reinforced with steel fibres. The aim of the tests was to determine the fracture characteristics of two types of interfaces. Beams of dimensions: 150 × 150 × 550 mm with a notch of 50 mm high were investigated. For each concrete, the cubic compressive ( f cm,cube ), flexural ( f ctm, f lex ), and splitting ( f ctm,split ) tensile strength were determined as well as fracture parameters like the fracture energy (G f ) and the fracture toughness. The most important results of these laboratory tests, important for the present studies, are summarized in Table 1. The bond effectiveness coefficient α int for the interface can be estimated at 0.46. The values of fracture energies seemed to be quite high-their values were greater than the ones predicted by code formulas [35,36]. The obtained values of the fracture energy for homogenous specimens were also greater than the one reported by Shah and Chandra Kishen [37]. They obtained G f = 157 ÷ 172 N/m for the concrete of similar strength as the low-strength concrete ( f cm = 34 MPa) and using a very similar testing procedure. On the other hand, values G f for the interfaces were close to each other (≈70 ÷ 80 N/m in reference [37]).  [26]. n/a-not applicable, n/t-not tested, n/g-not given.

Own Splitting and Shear-Slant Tests
A pilot study, which facilitates the determination strength of the interface, was performed at the Faculty of Civil Engineering WUT. The research program consisted of slant-shear and splitting tests. The slant-shear specimens were of a rectangular shape with dimensions 100 × 100 × 500 mm, and the shear plane was placed at 30 degrees with the vertical direction. The splitting tests were performed using cubic specimens of 150 mm edge. The assumed class for the stronger concrete was C30/37 (concrete A) and for the stronger one was C55/67 (concrete B). The mix composition was summarized in Table 2. During each casting, six additional 150 mm cubic specimens were prepared in order to determine the cubic compressive (three specimens) and splitting tensile strength (three specimens). The cubic compressive and tensile strengths were determined in accordance with their standards, respectively, [38,39]. Firstly, concrete B was cast (substrate), then, after 7 days, concrete A was poured into the mould (overlay). Before casting the weaker concrete, the contact surface of the hardened half of the specimen was subjected to chipping with a manual hammer. The specimens were cured in water up to the strength tests in order to reduce the influence of shrinkage. The strength tests were performed after 28 days since the second casting. The results of the tests are summarised in Tables 3 and 4.
The mean values, the standard deviation (STD), and the coefficient of variation (CoV) were reported. The specimens after splitting and shear-slant tests are shown in Figure 3. In case of slant-shear, the pure adhesive failure mode was obtained for all specimens, whereas the specimens in splitting the test presented a mixed failure mode (adhesive-cohesive). The bond efficiency coefficient was quite high (0.88) due to the early age of the substrate concrete and the surface preparation technique. The coefficient of variation for the tests of the interface was much higher than for homogenous specimens, which is consistent with other reported research data [16,25].
The following well-known formula for determining the splitting tensile strength was used: where: F ult -the ultimate force, L = 150 mm-the length of load line, and d = 150 mmthe length of sample destruction cross-section. The relations between the ultimate load and strength characteristics of an interface (σ ult -the ultimate compressive normal traction, τ ultimate -the ultimate shear traction) are as follows [25]: where: α = 30 • -the interface angle with the vertical, and a = 100 mm-the length of the cross-section edge.

Splitting and Shear Slant Tests by Santos and Julio (2009, 2011)
Santos and Julio performed the wide laboratory test programme, which aimed at assessing the shear strength between concrete layers [21,25]. They analysed the influence of factors such as the surface preparation, shrinkage, the difference of ages, and concrete stiffness. For the purpose of the present studies, specimens from series L28 and substrate surface-treated with wire-brushing (WB) were chosen. In this series, the added concrete was cast after 28 days. The same two types of tests were performed like in the pilot studies-the splitting and slant-shear tests. To determine the compressive (three cubes) and splitting tensile (five cubes) strength, standard cubes of edge length equal to 150 mm and five specimens with dimensions 150 × 150 × 450 mm were used during slant-shear tests. The splitting test was modified by using a flat surface in order to load two parts of composite specimens more evenly. The most important outcomes of the laboratory tests are summarised in Tables 5 and 6.

Numerical Modelling Strategy
The interface was modelled using cohesive elements, whose symbol is COH3D8 according to programme documentation [15]-see Figure 4. Such elements have eight nodes, so they use a linear approximation of the displacement field inside the element. They can be used with continuum-based constitutive relations or with a constitutive model formulated in terms of a traction-separation function. In a 3D domain in the elastic regime (before damage initiation), the traction-separation law neglecting normal-shear states coupling has the following form: where: t-the traction stress vector; indexes n, s, and t mean normal, first, and second tangential direction, respectively, e.g., t n -traction normal to the interface; K nn , K ss , and K tt -the stiffness of interface material in three directions, n = δ n T 0 , s = δ s T 0 , and t = δ t T 0 , respectively-interface nominal strains, T 0 -initial thickness of the interface, δ n , δ s , and δ t -separations (displacement jumps at an interface) in three orthogonal directions. It is worth mentioning that the traction stress vector is related to the Cauchy stress tensor (σ) in the following manner t = σ n, where n is normal to the interface plane [40,41].
Equation (4) can be rewritten as: where: K nn = K nn T 0 , K ss = K ss T 0 , and K tt = K tt T 0 -which are stiffness factors of the interface. In such a form, the elastic part of constitutive model for interface is given in references [13,40]. The stiffness coefficients for concrete-concrete interfaces have no clear physical interpretation since the interface before degradation can be assumed as a rigid one. Some authors try to prescribe a T 0 physical interpretation, e.g., the interface zone was assumed as a 100 µm thickness in reference [33]. However, assuming very small interface thickness, resulting in large values of its stiffness, can have an adverse influence on the convergence of the incremental-iterative procedure due to the large values of the unbalanced force vector. On the other hand, the introduction of too-small values can result in the wrong prediction of stress concentrations, which occur in the vicinity of two material connections. Therefore, the recommendations from the ATENA manual seem to be reasonable [13]. It recommend there to assume the following stiffness for the cohesive elements for the concrete-concrete interface: where: a-the biggest dimension of the connected parts and E and G are Young's and Kirchhoff's modules of the weaker concrete, respectively. The inelastic part of the traction-separation law implemented in the Abaqus code was formulated in the framework of the damage theory. The onset to the damage state is indicated by the damage initiation criterion. Among the few available ones, the quadratic nominal stress criterion was selected. It can be represented as: where: t max n -maximum normal traction stress (i.e., tensile strength, associated with mode I in the nomenclature of fracture mechanics); t max s and t max t -maximal shear tractions in directions n and s (i.e., shear strength in pure shear), respectively. . . . -the Macaulay brackets, which means the following operation: Therefore, damage factors are not activated in a pure compression state, and allowable shear stresses do not depend on normal stress-see Figure 5. Such behaviour does not correspond to the concrete-concrete interface, since compressive normal stress causes an increase in the shear strength. Carol's formula [42] is believed to describe the failure envelope for the concrete-concrete interface correctly [18,40,43]: where: τ = t 2 s + t 2 t -the resultant shear traction; f t -the tensile strength of interface; and c-the cohesion. The standard Abaqus damage traction-separation interface model was enhanced by a user subroutine USDFLD to make the failure envelope pressure-dependent: where: , f ctm and f sh -the tensile and shear strength, respectively. The USDFLD subroutine allowed us to modify the parameters of models implemented in the Abaqus code by default at the beginning of each increment of the Newton scheme. Consequently, it is important to set a small value of increment load size when one uses such an approach. On the other hand, the problem is highly non-linear, so small increments are necessary for the convergence of the incremental-iterative procedure. The proposed algorithm reads the normal traction and calculates the admissible shear stress according to Carol's formula. All the discussed criteria are shown in Figure 5.
The last issue connected with a traction-separation interface model is its post-cracking behaviour. After reaching the damage-initiation criterion, tractions are calculated taking into account the damage factor D: The above equation shows that in the case of interface compression, there is no reduction in normal stiffness. The evolution of the damage parameter under a combination of normal and shear deformation is described using an effective displacement: In the present modelling strategy, exponential damage evolution was used, in which damage variable D is described by the relation: where: δ max m -the maximum effective displacement during the loading history; δ init mthe effective displacement at the damage initiation; δ f ail m -the effective displacement at complete failure; and α-the parameter which controls the rate of softening (in the present studies, the typical assumed value for concrete was equal to 7.0 [44]). The model behaviour is shown in Figure 6.
The post-peak behaviour of the concrete-to-concrete interface differed significantly for tension and shear states. In the present study, the effective displacement at failure was related to the fracture energy for the mode I type of fracture (G f ,I ) or for mode II (G f ,I I ). The energy dissipated in any mode (G f ) can be calculated according to: where: t ult -the ultimate stress in the given mode. The first term can be neglected due to its very small value in comparison to the second one. After calculating the integral and some mathematical manipulations, the formula for the effective displacement of failure can be obtained: The fracture energy in mode I is similar to the energy in mixed mode I-II conditions [45], whereas in a pure mode II, it is much larger-about 20 times more [46]. The fracture energy, when the interface is under simultaneous compression and shear, is even greater, since, after the crack formation, the interface is capable of transmitting shear stresses due to the shear-friction mechanism [42,47,48]. Tests for a concrete-to-concrete interface, in which the fracture energy was determined in addition to its strength, are scarce, so the assumption of the correct value of this parameter is not an easy task. The following relation can be assumed between the fracture energy of the interface (G f ,int ), the bond efficiency coefficient (α int ), and the fracture energy of rgw weaker concrete (G f ,weak ): where: β int -the coefficient that can be related to α int . These relations obtained by three research teams ( [26,34,37]) are summarised in Figure 7. After analysing the data, it can be concluded that the value of coefficient β should be taken from 1.2 to 2.0. In case of a lack of experimental results, the fracture energy of weaker concrete in modes I and II can be estimated with the following formulas [36,46]: where: d max -the maximum aggregate size (in (mm)), and f cm,cube -the cubic compressive strength (in (MPa)). The whole calibration strategy for a concrete-to-concrete interface was summarised in Table 7. The viscous regularisation was used to obtain a stable model response in the post-peak regime. It is worth mentioning that in the CAE module of the Abaqus code, this parameter has to be introduced via the "Assign Element Type" window, which is not in the "Property" module. µ Viscosity parameter 0.0001 The concrete region was modelled with C3D8 continuum elements with selective integration [15]. The concrete damage plasticity (CDP) was chosen as a constitutive model for concrete. The version without the scalar damage parameter was assumed since only monotonic loads were analysed. The theoretical background of this model was presented in references [50,51]. The model has gained significant popularity and was discussed in detail in many references, e.g., [52][53][54][55], so its thorough description could be omitted in this study.
In general, the model is formulated in the framework of the theory of plasticity with an additive strain tensor decomposition: where: and pl -total and plastic strain tensors, respectively, and D el -the elasticity tensor for isotropic material. The CDP model uses a yield criterion proposed by Lubliner [50] and isotropic strain hardening/softening. The flow is governed by a non-associated flow rule. As the plastic potential, the smoothed Drucker-Prager cone is assumed. The hard-ening/softening in tension is controlled by equivalent plastic strain in tension (PEEQT according to the manual [15]), whereas, in compression, another equivalent plastic strain is introduced (PEEQ according to [15]). Distribution of these scalars can be used to analyse which regions of the model are cracked (PEEQT) and which regions are crushed under compression (PEEQ) during the loading history of the analysed specimen. The parameters that have to be entered into the programme are summarized in Table 8.
In all presented case studies, the shape of the loading device was modelled precisely. In order to reduce the computational effort, they were discretized with R3D4 rigid elements. Between the surfaces of the loading device and the specimens, the contact with Coulomb friction was assumed with the friction coefficient equal to 0.6. Due to convergence issues, the simulations were displacement-driven, i.e., as the boundary conditions, displacements of loading devices were assumed.
The non-linear problem was solved using an incremental-iterative Newton-Raphson algorithm. Two default convergence criteria were used, which are precisely described in reference [3]. The tolerance for unbalanced forces was assumed as equal to 2%, whereas the tolerance for displacement correction vector was set to 1% (default values). µ Viscosity parameter 0.0001

CS1-Tests with One Cohesive Element
The first example concerns a numerical model consisting of one C3D8 continuum and one COH3D8 interface element-see Figure 8. Mechanical parameters adopted for constitutive models can be found in Table 9. Results for three loading paths are presented: • displacement-driven tension (pure mode I), i.e., δ n > 0, δ s = δ t = 0, • displacement-driven pure shear (pure mode II), i.e., δ n = δ t = 0, δ s > 0, • at the first step, the pressure that induces compressive normal traction; at second-step displacement-driven shear, i.e., δ n = δ t = 0, δ s > 0.
The results of all performed analyses are shown in Figure 8b in one plot to highlight the differences in the dissipated energy in mode I and mode II. Moreover, it is clearly visible that the model predicts an increase in the shear strength due to a normal compressive traction. Some discrepancy between the model prediction and the real behaviour of a concrete interface is the softening branch shape in compressive states. Usually, it is assumed that the interface has residual strength in the case of a shear under compression due to shear forces [40,42,57], whereas the model predicts that shear traction stresses tend to be zero. On the other hand, the additional amount of energy is dissipated in this case, which may be related to shear forces. Summarising, the model behaviour is consistent with the assumptions from Section 3.

CS2-Simulation of Three-Point Bending of a Notched Beam
The next example is a simulation of a three-point bending test of a beam with a notch described in Section 2.1. Due to the symmetry of the specimen, half of the beam was modelled. The qeometry of the model with assumed boundary conditions (symmetry condition on the symmetry plane and simple supports) is shown in Figure 9. Mesh dependency studies were performed assuming three mesh densities in the vicinity of the notch and the interface: 10, 5, and 2 mm-see Figure 10.  The results of the performed analyses are shown in Figure 11 as force-crack-mouthopening-displacement (CMOD) curves. They were compared with the experimental outcomes reported in reference [26]. The best agreement was obtained for the densest mesh (2 mm). The little discrepancies may be caused by the non-homogenous distribution of concrete strength. In Figure 12, the distribution of normal traction in the interface is shown for three load levels (0.52 F ult , 0.96 F ult , and F ult ) as well as a deformed mesh.

CS3-Simulation of Own Tests
The following case study concerns the simulation of our own laboratory tests: splitting and slant-shear of a concrete prism. A detailed analysis of the one-material splitting test, sometimes called "a Brazilian test," can also be found in reference [58]. The models and meshes adopted in the analyses are shown in Figures 13 and 14. In the case of the splitting test, the double symmetry of the task was used-the quarter of the specimen was considered. As in the previous CS, three mesh densities next to the interface were studied: 10, 5, and 2 mm. The test stand was modelled precisely, e.g., the exact geometry of a loading device with a pad was reproduced. The ultimate forces obtained from laboratory tests and numerical analyses are shown in Figure 15. The model for the splitting test slightly underestimated the ultimate force, whereas the model for the slant-shear test predicted value that was slightly too large. However, the differences were not too great (4-8%). A deformed mesh and a map of equivalent plastic strain in tension (PEEQT according to the programme documentation [15]) for the splitting probe are shown in Figure 16. The plastic strains occurred in the weaker concrete, so the stresses during the analysis exceeded its tensile strength. This observation explains the increased possibility of mixed failure mode (adhesive-cohesive) in the case of the interface with a large value of α int loaded in accordance with the standard procedure [39,59]. This is caused by problems with the uniform loading of both parts of the bi-material specimen.
A deformed mesh and map of equivalent plastic strain in compression (PEEQ according to the programme documentation [15]) for the slant-shear test are shown in Figure 17. The failure mode was reproduced correctly. Concentrations of plastic strains near the sharp edge explain the characteristic chips that appeared after the strength test of such specimens.  Figure 18 summarises the result of the parametric studies performed using the slantshear test model. They were performed in order to asses the influence of the viscous parameter µ, which has no clear physical meaning and which was introduced to stabilise the incremental-iterative process. The fracture energy of an interface (G f ,I and G f ,I I were equally increased/decreased) was estimated according to Section 3. The parametric studies revealed that parameter µ has a very small impact on the results. Conversely, the ultimate force prediction strongly depends on the fracture energy, which proves that the issue of its proper estimation is of great importance.

CS4-Simulation of Santos and Julio Tests
The last CS concerns the simulation of Santos and Julio [21,25] tests described in Section 2.3. Since our own tests followed the same testing programme, the boundary conditions applied in the analysis were almost the same except for the method of splitting the specimen loading. The specimen was loaded with a plain, narrow metal surface, not with a cylindrical-shaped device-see Figure 19a. The mesh of size 5 mm was used because of a small change in the ultimate force prediction in comparison to a 2 mm mesh in the previous CS. The mesh for a slant-shear test is shown in Figure 19b because of a different specimen size (CS3-100 mm × 100 mm × 500 mm; CS4-150 mm × 150 mm × 450 mm). A comparison of results obtained with the numerical analysis and from experimental tests is shown in Figure 20. The agreement between them was very good (with differences not greater than 2%), which proves the reliability of the formulated modelling strategy. In Figure 21, the evolution of the scalar degradation parameter D for an interface is presented as maps for three load levels: 0.6 F ult (initiation of damage), 0.94 F ult , and F ult . Traction distributions along the selected cross-section (in the middle of the specimen) are shown in Figure 22. In the elastic stage (0.5 F ult ), stress concentrations next to the specimen edges, caused by its inhomogeneity, can be observed. For the ultimate force (F ult ), the delamination of the specimen is progressing (decrease in tractions at the left edge). In the post-peak regime, the rest of the concrete contact is additionally loaded (increase in tractions), so the proper definition of the strength envelope is crucial also in this stage. It is worth mentioning that in real structures the progressive delamination of two concrete interfaces rarely leads to reaching the load capacity of the whole structural system.

Discussion
The strategy for adjusting an interface model available by default in the Abaqus code to the behaviour of concrete-to-concrete interface was proposed in Section 3. It was based on cohesive elements with the damage-elastic traction-separation constitutive model, which should be used to model initially rigid and thin connections between layers made of different concretes [13,15,40]. The standard model was enhanced using the simple USDFLD user procedure, which enabled us to introduce the strength envelope, which is dependent on a normal traction value. The model took into account the different fracture energy values for modes of fractures I and II. The proposed approach was verified and validated with four case studies concerning: one element test, simulation of three-point bending of the bi-material notched beam [26], simulation of two series of splitting, and slant-shear tests-our own pilot research and the research made by Santos and Julio [21,25]. The results of the performed analyses are summarised in Table 10. The ratio of the ultimate force predicted by the FEM model and obtained as an outcome of a laboratory test was calculated for CS2-CS4. Its mean value was 1.01, and the CoV was 5%, which proves the the accuracy of the proposed modelling strategy.
Besides, the fully non-linear simulation of the most popular laboratory tests aimed at determining the interface strength gave a better insight in the specimen's behaviour under loading. The simulation for the three-point bending test of the notched beam showed the progressive delamination in the interface, cf. Figure 12. Such a phenomenon can be examined using modern optical measurement techniques [27]. On the contrary, such techniques fail in monitoring fracture process for other analysed tests (the splitting and slant-shear tests) due to the brittle failure mode of such specimens [23]. Therefore, the numerical methods are appropriate tools to understand their behaviour in depth. In the present research, the characteristic features of their failure modes were reproduced, such as the mixed failure mode for cubic specimens in the splitting test-see Figure 16-as well as the formation of chipping regions in the slant-shear test-see Figure 17. Moreover, the traction distributions along the interface in different stages of the tests were shownsee Figure 22.
The limitation of the presented approach is the inability of the model to cover the residual strength envelope caused by the shear-friction phenomenon and the irreversible (plastic) slip since the default traction-separation material model was formulated in the framework of damage-elasticity without the plastic part. These issues can be overcome by the implementation of the traction-separation law within the UMAT or UEL user procedure [57].

Conclusions
The article presented the complete strategy for calibration of the concrete-to-concrete interface model in the Abaqus code. The default traction-separation implemented in this system does not take into account the strength envelope for such interfaces. This issue was overcome with the USDFLD user subroutine coded in FORTRAN, which enabled us to relate a shear strength with tractions normal to the interface. Different values of fracture energy in fracture modes I and II were also included. The results presented in the article can be used in the analyses of complex composite structures made of concrete layers cast at different times as well as in the assessment of various methods of repairing deteriorated structures. The presented modelling approach can also be useful in the simulation of masonry structures-the interface model can be easily adjusted to the behaviour of mortar [41].
Furthermore, very detailed simulation of the most popular strength tests' methods for concrete-to-concrete interfaces was presented in four case studies. The following observations can be made on their basis:

•
The method of loading specimens, which is fully consistent with the standards for splitting tests [39], can increase the possibility of obtaining cohesive or mixed failure modes (see Figure 16), so the modification of the test procedure introduced by Santos and Julio [25] seems to be reasonable. • In the case of small specimens, the mesh density of 5 mm was sufficient to obtain satisfactory accuracy of the FEA results. • The numerical model was able to cover the chipping of sharp edges for slant-shear specimens (see Figure 17). • The value of fracture energy has a noticeable influence on the ultimate force predicted by the model (see Figure 18), so this issue should be experimentally studied in more depth due to the small amount of the experimental data. • The traction stress distribution is not homogenous along the interface during the whole loading history in the case of slant-shear specimens made of concrete, which are of different classes (see Figure 22 and reference [21]). Consequently, the interface strength characteristics determined according to Equations (3) should be corrected in order to take into account stress concentrations present in the prism specimens.
Further studies will be concentrated on the determination of the fracture energy for different concrete-to-concrete interfaces (different concrete classes and surface preparation) and its dependency on the bond efficiency coefficient. Additionally, an attempt will be made to implement more of the sophisticated interface models-examples of such models can be found in references [40,42,57].

Data Availability Statement:
The data presented in this study are available on request from the authors.