XFEM-Based Multiscale Simulation on Monotonic and Hysteretic Behavior of Reinforced-Concrete Columns

Abstract

The extended finite element method (XFEM) is efficient in simulating crack initiation and its evolution process for reinforced-concrete (RC) structures due to its ability to solve fracture problems. Moreover, the multiscale numerical simulation helps understand global and local failure behavior of RC structures simultaneously. In this study, the XFEM-based multiscale modeling approach was proposed to investigate the monotonic and hysteretic performance of RC columns. Firstly, two-scale models composed of fiber beam elements and XFEM-based solid elements with homogeneous material assumptions were established using compiled material subroutines for fiber beam elements. Secondly, the accuracy of XFEM-based two-scale analysis in predicting the hysteretic behavior of tested RC columns was verified by comparing the crack morphology and load-displacement curve obtained from tested specimens under different axial compression ratios (ACRs) and two-scale models using the concrete damaged plasticity (CDP) model. Thirdly, multiscale models of RC columns were constructed with fiber beam elements, XFEM-based solid elements and mesoscopic concrete models composed of mortar, interfacial transition zone (ITZ) and aggregates with different geometric shapes and distribution patterns. Finally, the XFEM-based multiscale simulation was employed to investigate the influence of mesoscale structure variation of concrete on both global behavior and local failure patterns of RC columns subjected to monotonic loading. The simulation results of multiscale models established with CDP model and XFEM were comparatively discussed in depth. The XFEM-based multiscale simulation developed in this study provides an efficient modeling approach for investigating the stochastic nature of cracking behavior in RC columns.

1. Introduction

In order to clearly understand the fracture mechanism of reinforced-concrete (RC) structures, efficient numerical simulation on the whole failure process of RC structures from the initiation of cracks to the final failure at both local and global points of view is meaningful. For extending the simulation accuracy of a traditional finite element method (FEM) in simulating cracks of concrete materials and RC structures, various cracking models for concrete materials have been developed [1,2]. In particular, the concrete damaged plasticity (CDP) model has been extensively employed to simulate the damage initiation and evolution process in concrete [3,4], where the damage coefficient is adopted to reflect the degradation of material strength and stiffness of cracked concrete members. However, the concrete is a typical elastic-brittle material, and it changes from a continuum component to a discontinuous body once a crack initiates, which leads to degradation in material strength, structural stiffness and load-carrying capacity of RC members. It is difficult to truly reflect the crack initiation and propagation using damage coefficients.

In recent years, the extended FEM (XFEM) has been proposed to describe the discontinuous fields independent of element boundary due to its strong vitality in analyzing discontinuous problems. Utilizing XFEM, the cracking process of the simply supported beam was simulated with preset cracks [5]. Yin et al. [6] also introduced XFEM to the numerical analysis of concrete slabs and beams. In addition, the pinching mechanism of RC columns under cyclic loading has been systematically investigated as well [7]. Moreover, the strain fluctuation and stress redistribution of steel bars and concrete in RC structures subjected to cyclic loadings can be obtained [8]. The opening and closing phenomenon of each crack can be effectively simulated at macro-scale [7,9]. XFEM has become a popular method to simulate crack propagation problems in concrete [10]. However, it has been recognized that the nonuniformity of concrete at meso-level and the randomness of aggregate distribution in concrete are essential factors leading to the non-linear properties of concrete materials and the stochastic structural behavior at macro-level [11,12,13,14]. Until now, perfect simulation for predicting the local crack evolution and failure patterns of concrete are still challenging since the properties of concrete materials are mainly obtained from material tests at macro-level [1]. At meso-level, concrete is treated as composite material composed of aggregate, mortar, interfacial transition zone (ITZ), initial pores and defects. By means of XFEM, the formation and propagation processes of microfracture in concrete have been studied at meso-level [15], which laid a solid foundation for XFEM-based mesoscale simulation of concrete.

The two-scale model composed of beam elements and solid elements or shell elements has been proposed and extensively employed in FEM analysis on RC structures owing to its high computational efficiency [16]. Meanwhile, various kinds of two-scale modeling methods have been successfully applied in the damage and failure analysis of civil structures under earthquakes [16,17,18,19,20]. Moreover, the mesoscale modeling approach has been successfully applied in the numerical simulation of concrete structures subjected to various loading cases [21,22,23]. The feasibility of XFEM-based multiscale models composed of two-scale and mesoscale models needs to be further investigated to provide an efficient way to analyze the structural behavior of RC structures.

In this study, the XFEM-based multiscale numerical modeling methodology was developed to investigate both global and local failure behavior of RC columns. With the help of XFEM, two-scale models composed of fiber beam elements, solid elements and multiscale models composed of fiber beam elements, solid elements and mesoscale concrete constructed with circular, polygonal and elliptical aggregates with various distribution patterns, mortar, and ITZ were established to study the crack behavior of tested RC columns under monotonic and hysteretic loadings, respectively. Two material subroutines for fiber beam elements suitable for implicit analysis and explicit analysis of RC columns were developed. The accuracy of developed material subroutines was validated using experimental results and multiscale simulation results with a CDP model. The monotonic and hysteretic behaviors of RC columns under different axial compression ratios (ACRs) were numerically investigated. Finally, the XFEM-based multiscale models were employed to evaluate the influences of mesostructure variation of concrete on both global behavior and local failure patterns of RC columns. The advantages and disadvantages of XFEM and CDP-based multiscale modeling approaches were discussed comparatively.

2. Description of Tested Reinforced-Concrete (RC) Columns

2.1. Specimen Layout

The mechanical behavior of four RC column specimens subjected to monotonic and cyclic horizontal loadings was tested. The height of RC columns was 1150 mm, and the cross-sectional dimension of tested RC columns was 250 mm by 250 mm. The protective layer thickness was set to 25 mm. The diameters of stirrups and longitudinal steel bars were 6.0 mm and 12.0 mm, respectively. The distance between stirrups was 50 mm, and it was set to 30 mm in the vicinity of the loading point to avoid local failure at the top of the tested columns. The columns were connected with the stiff floor of the laboratory. There RC column specimens subjected to cyclic loading were labeled as C1 to C3, which were loaded with axial compression ratios (ACRs) of 0.08, 0.10 and 0.20, respectively. The tested specimen under monotonic loading with the ACR of 0.08 was numbered as C4. The geometry dimension and reinforcement arrangement of C1–C4 were identical, as shown in Figure 1.

2.2. Material Properties

The concrete grade was originally designed as C40 in the China design code. The concrete mix design is detailed in

Table 1. Concrete mix design (kg/m³).

Portland Cement Grade 42.5	Water	Medium Coarse Sand	Aggregates	Water Cement Ratio	Sand Ratio (%)
474	180	476	1320	0.38	26.5

. The material of steel reinforcement was HRB 335 (hot-rolled ribbed steel bars) and HPB 235 (hot-rolled plain steel bars). According to experimental observation from material strength tests, the compressive strength of concrete was set to 40.9 MPa for numerical simulation. Moreover, the tensile strengths for stirrups and longitudinal bars were 551.4 Mpa and 519.0 Mpa, respectively. The yielding strengths were set to 378.9 Mpa and 371.5 Mpa.

2.3. Test Setup and Loading Protocol

The vertical axial force was applied on the top of RC columns by a hydraulic jack, and the horizontal loading was applied by an actuator at the height of 1.0 m. The horizontal displacement at the same height of the loading point was measured with linear variable differential transformers (LVDTs). The vertical load was measured with a pressure gauge between the rolling support and hydraulic jack, as shown in Figure 2. In the experiment study, the vertical load was applied by load control, and then kept constant during the whole loading process. The horizontal load was directly applied under displacement control. The loading method in the numerical simulation was consistent with that of tests.

In this study, the measurement devices included the strain gauges, displacement meters, dynamometer and the data acquisition system. For concrete, the resistance value of strain gauges was 120 Ω ± 0.2%, with the sensitivity coefficient of 2.032 ± 0.26%. For steel bars, the resistance value of strain gauges was 119.8 ± 0.1 Ω, with the sensitivity coefficient of 2.08 ± 1%. The relative errors of LVDT and dynamometer in MTS loading system and hydraulic jack were both ±0.05%.

3. Modeling Procedures of Extended Finite Element Method (XFEM)-Based Multiscale Analysis

3.1. Theoretical Foundation of XFEM

The XFEM is based on the concept of partition of unity method (PUM) [9], which is an extension of conventional finite element methods. The singularity around the crack tip is described by a near-tip asymptotic function. A discontinuous Heaviside function is used to express the displacement jump at the fracture surface. The enrichment function is composed of near-tip asymptotic functions and discontinuous function. After integrating the partition of unity enrichment, the displacement vector function $u$ can be approximately expressed as Equation (1) [24].

$$u={\displaystyle \sum }_{I=1}^N{N}_I\left(x\right)[u_I+H\left(x\right)a_I+{\displaystyle \sum }_{a=1}^4{F}_{\alpha }\left(x\right)b_I^a$$

(1)

where $N_I\left(x\right)$ is the shape function of node $I$ and $u_I$ stands for the displacement vector for the continuous part without cracks passing through the crack tip. $a_I$ and $H\left(x\right)$ represent the freedom vector of an enriched degree and the Heaviside function. $b_I^a$ is the degree of freedom (DOF) related to elastic asymptotic crack-tip functions $F_{\alpha }\left(x\right)$. $u_I$ is suitable for all nodes. $H\left(x\right)a_I$ is only applicable for nodes where the element is passed through by cracks. ${\sum }_{a=1}^4{F}_{\alpha }\left(x\right)b_I^a$ is suitable to nodes where the corresponding element contains the crack-tip. When the element is not cut by cracks and does not contain crack-tip, Equation (1) degenerates into the displacement formula of conventional FEM. The normal and tangential coordinates for a smooth crack are illustrated in Figure 3.

The generalized Heaviside function is employed to describe the displacement discontinuity due to cracks. The mathematical expression of $H\left(x\right)$ can be described as Equation (2).

$$H\left(x\right)=\{\begin{array}{ll}+1& \left(x-x^{*}\right)·n\ge 0\\ -1& others\end{array}$$

(2)

where $x$ stands for the sample point for investigation, $x^{*}$ refers to the point on crack which is closest to $x$. $n$ represents the unit outward normal to crack at point $x^{*}$.

The asymptotic crack tip function $F_{\alpha }\left(x\right)$ can be described as Equation (3)

$$\left[F_{\alpha }\left(x\right)\right]=\sqrt{r}\left[sin\frac{\theta }{2},cos\frac{\theta }{2},sin \theta cos \frac{\theta }{2},cos\theta sin\frac{\theta }{2}\right]$$

(3)

where $r$ and $\theta$ represent the polar coordinates. At crack tip, $\theta$ is equal to zero.

3.2. Multiscale Modeling Methodology

In this study, the multiscale modeling methodology for RC columns is proposed. The flowchart of modeling methodology for multiscale simulation is presented in Figure 4. In general, the fiber beam element model is efficient in analyzing the global structural performance of RC columns. However, the detailed damage evolution process cannot be obtained. In order to describe the damage initiation and development at the foot of RC columns, solid elements are usually preferred in numerical analysis. Since the failure of RC columns usually occurs at the bottom of tested specimens, two-scale models are employed to investigate the global behavior and local failure pattern of tested specimens. For improving the computational efficiency, fiber beam elements are employed to model the upper part of RC columns. It can be seen that the two-scale model is the combination of the fiber element model and solid element model [16]. Multiscale models are established by modeling the lower part of two-scale models with mesoscale concrete.

As shown in Figure 5, the multiscale concrete models are composed of fiber beam elements, solid elements, and mesoscale concrete with different aggregate distribution patterns and geometric shapes. For comparison, the CDP model is adopted for the material definition of mortar and the ITZ, as presented in Figure 5a. Moreover, the aggregates are assumed as elastic components. In Figure 5b, the solid element segment and mesoscale concrete segment are constructed using XFEM. Since the ITZ layers have been explicitly modeled, the number of elements increases significantly. Without loss of generality, five sets of random aggregate samples are generated to investigate the effect of mesostructure variation of concrete on local and global behaviors of RC columns.

Moreover, the separation of fiber element part and solid element part in two-scale models should be determined reasonably to balance computation efficiency and accuracy. For fiber beam elements, the constitutive models for steel bars and concrete are compiled using user material subroutines (Umat and Vumat) in ABAQUS, as detailed in Section 3.3. The height of the RC column modeled with solid elements $l_{solid}$ can be determined according to the plastic segment length $l_p$ calculated using Equation (4) [25].

$$l_p=0.08L+0.022d_b{f}_y$$

(4)

where $d_b$ and $f_y$ stand for the diameter and yield strength of longitudinal reinforcement, respectively.

Furthermore, the multiscale models composed of fiber beam element, solid element and mesoscale concrete are constructed to consider the mesostructure variation of concrete. According to the practical crack distribution and failure pattern of tested RC columns, the height of the mesoscale concrete segment is set to $0.2 L$.

In this section, the behaviors of RC columns under monotonic and cyclic loadings, which are simulated using the fiber beam element model, solid element model and two-scale model composed of fiber beam elements and solid elements, are compared to verify the feasibility of the developed two-scale modeling approach. As illustrated in Figure 5a, the CDP-based two-scale and multiscale model are also established to further validate the accuracy of developed material subroutines for the fiber beam element and XFEM-based multiscale modeling approach. The parameters of the CDP model are shown in

Table 2. Material parameters.

Material	Young’s Modulus (GPa)	Poisson’s Ratio	Compressive Strength (MPa)	Density (kg/m3)	Fracture Energy (N/m)	Maximum Principal Stress (MPa)
Aggregates	50.0	0.16	/	2733	288.0	6.0
Mortar	24.7	0.20	27.2	2400	168.0	3.5
ITZ	21.2	0.20	19.3	2400	115.2	2.4
Concrete	31.3	0.20	40.9 *	2400	216.0	4.5
Stirrup	200.0	0.30	371.5 * (yielding stress)	7800	---	---
Longitudinal steel bar	200.0	0.30	379.6 * (yielding stress)	7800	---	---

Note: * denotes the material parameters direly measured from material strength tests.

. The damage coefficients of the CDP model are calculated according to Equations (5) and (6) proposed by Birtel and Mark [26]:

$$d_c=1-\frac{{\sigma }_c/E_c}{{\epsilon }_c^{pl}\left(1/b_c-1\right)+{\sigma }_c/E_c}$$

(5)

$$d_t=1-\frac{{\sigma }_t/E_c}{{\epsilon }_t^{pl}\left(1/b_t-1\right)+{\sigma }_t/E_c}$$

(6)

where $E_c$ is the initial modulus, ${\sigma }_c$ and ${\sigma }_t$ stand for the compressive and tensile strengths of concrete. ${\epsilon }_c^{pl}$ and ${\epsilon }_t^{pl}$ are equivalent compressive plastic strain and equivalent tensile plastic strain.

The XFEM-based multiscale simulation proposed in this study is exhibited in Figure 5b. By contrast with the two-scale modeling approach with conventional solid elements, the two-scale model herein is composed of XFEM-based solid elements and fiber elements. The multiscale models of RC columns are established with fiber beam elements, XFEM-based solid elements and XFEM-based mesoscale concrete constructed with mortar, the ITZ and randomly distributed aggregates. The material properties of the steel reinforcement, the XFEM parameters of homogeneous concrete, the ITZ, mortar and the aggregates are detailed in Table 2.

For reducing computational complexity and avoiding severe convergence problems in 3D XFEM, the numerical models in Figure 5b are simplified as 2D cases. Similar to CDP-based two-scale and multiscale models, the interface between fiber elements and plane stress elements is connected utilizing a kinematic coupling method. The steel reinforcements are simulated with the truss element T2D2, and the four-node plane stress element CPS4R with reduced integration is adopted to model homogeneous concrete, mortar, ITZ and aggregates. In order to improve the calculation accuracy, the bond-slip between steel reinforcements and concrete is considered in XFEM models. The element nodes of steel bars and concrete are connected with spring connector elements. The bond-slip effect is neglected at the top half of RC columns.

It has been reported that the mesostructure variation of concrete leads to fluctuation of material strength and Young’s modulus of concrete subjected to axial tension and compression [3]. In order to investigate the influence of mesoscale structure of concrete on the cracking pattern of tested RC columns, the multiscale numerical models are constructed utilizing a random aggregate method (RAM). The modeling procedure of the RAM has been extensively reported [3,4,27], which is excluded in this section. Besides, the randomly-distributed aggregates are simulated using polygonal, circular and elliptical aggregates generated by the random aggregate generation and delivering program compiled in this study. As illustrated in Figure 5, three polygonal samples are designed to investigate the influence of distribution variation of aggregates. Meanwhile, elliptical and circular aggregate samples are considered to discuss the effect of aggregate geometry variation on the cracking evolution process and global structural response of tested RC columns. It is challenging to find a unified criterion for the definition of ITZ. The actual thickness of the ITZ is about 20–50 μm, and it is impossible to establish elements in numerical models with such a tiny thickness [28]. Taking 1.0 mm as the thickness of the ITZ is accepted in most of the literature related to mesoscale modeling of concrete [28]. Therefore, the ITZs in mesoscale modeling on RC components are simulated using one-layer elements with the thickness of 1.0 mm, as illustrated in Figure 5.

3.3. Material Subroutines for Fiber Beam Elements and Their Verification

A fiber beam element has been widely used in the numerical analysis for RC columns. As shown in Figure 5, the fiber beam element is employed for simulating the concrete and steel bars in the top half part of the two-scale model. The XFEM module has been implemented in ABAQUS. Moreover, the CDP constitutive law has been integrated into ABAQUS, which can be applied for 3D solid elements as well as 2D planar elements. For beam elements, there are no uniaxial hysteretic constitutive models available in ABAQUS. Therefore, the uniaxial hysteretic constitutive models for concrete and steel bars need to be developed for fiber beam elements employed in two-scale and multiscale models of RC column specimens in this study. The user-defined material subroutines have to be developed according to the interface requirements of material subroutine in ABAQUS. By referring to the hysteresis rule of Concrete 02 in Opensees [29], material subroutines for concrete are compiled in this study, as presented in Figure 6a. The hysteretic criterion for loading and unloading is defined according to the study performed by Yassin [29], which considers the stiffness degradation and hysteretic energy dissipation. To simulate the hysteretic behavior of steel bars, the bilinear model [29] is adopted as the constitutive law of steel reinforcement. In this model, the reverse loading curve directs to the point representing the maximum deformation during the loading history, as shown in Figure 6b. The Bauschinger effect is also considered in the compiled material subroutines.

The verification of material subroutines and simulation results with different modeling techniques is summarized in Figure 7. It can be seen that the numerical results of fiber beam models using Vumat and Umat are almost identical, as presented in Figure 7a,b. Figure 7c indicates that the simulated reaction force-displacement curves of the fiber beam model and two-scale model subjected to cyclic loading match well with that of solid beam element model and experimental result, which fully validates the accuracy of Umat and Vumat compiled in this study.

3.4. Material Properties of XFEM-Based Solid Elements

In the lower part of the two-scale model, the classic bilinear elastic-plastic model is adopted as the constitutive law of steel bars connected with solid elements. The elastic modulus, density and yield stress of the steel bars are detailed in Table 2. In order to consider the Bauschinger effect of steel under cyclic loading, kinematic hardening is adopted for material definition. The maximum principal stress criterion is used for homogeneous concrete, mortar, aggregate and the ITZ as the cracking initiation criterion of XFEM. The linear softening and B-K criterion are employed to define the damage evolution [1,7]. Since the compressive crush of concrete cannot be simulated using XFEM, the compressive stress–strain curve is defined as per the CDP model to consider the softening of concrete under compression [1,7]. The material parameters of both concrete and steel reinforcements involved in this study are shown in Table 2.

In conventional numerical simulation on RC structures, the concrete is usually assumed to be a homogeneous material without considering its mesostructures. It is relatively easier to define the material parameters for homogeneous concrete. However, the concrete at meso-level has to be modeled as composite component consisting of aggregates, ITZ and mortar, which makes it challenging to directly obtain all material properties for each phase from material tests. Similar issues also exist in the mesoscale modeling of concrete using the XFEM method. The determination and calibration of fracture energy and the maximum principal stress of ITZ are also critial issues. Therefore, the material parameters which cannot be easily obtained from material strength tests were determined as per the definition rules in existing mesoscale simulation on RC structures.

3.5. Boundary and Convergence Setting

As shown in Figure 5, all of the DOFs at the bottom of the RC column are constrained and the vertical loads are applied at the top of tested specimens. The monotonic and cyclic loadings with displacement control are adopted both in experimental study and numerical analysis.

3.6. Bond-Slip Effect and Spring Connector Elements

As present in Figure 5, the bond-slip effect between concrete and steel reinforcement is considered. The theoretical model of the bond stress–slip relationship reported in the Europe concrete code is employed [30], where $s_1=s_2=0.6 \mathrm{mm}$, $s_3=1.0 \mathrm{mm}$, $\alpha =0.4$, ${\tau }_{max}=2.0\sqrt{f_{ck}}$ and ${\tau }_f=0.15{\tau }_{max}$, as presented in Figure 8. However, a large number of non-linear connection springs can easily lead to severe convergence problems. In order to simplify the FEM models, the bond-slip between stirrups and concrete is neglected in this study. The spring stiffness of the connector element in the longitudinal direction is defined according to the analytical bond stress-slip relationship presented in Figure 8. The connection between concrete and longitudinal steel bars in their lateral cross-section plane is simulated using rigid spring connectors.

In general, in order to improve calculation efficiency, a reliable bonding condition between the steel reinforcements and the concrete is usually assumed in most FEM analysis. In this section, the influence of bond-slip on the hysteretic behavior of RC columns is investigated. Figure 9 shows the comparison of the hysteretic behavior of specimen C3 with and without considering the bond-slip in a solid element model and two-scale model. As shown in Figure 9, the pinching effect in hysteresis loops becomes more obvious and the structural stiffness of RC columns is obviously reduced while considering the bond-slip effect. When steel bars are directly embedded into concrete, the peak load increases by 8.6%. Crack distribution patterns corresponding to a solid element model and two-scale model are shown in Figure 9. The crack patterns present limited differences in either the solid element model or two-scale model whether bond-slip is considered or not. However, the specific location of the maximum fracture is slightly different. When considering bond-slip behavior, the maximum fracture appears at the column foot. The position of the maximum crack ascends and is mainly located around 8–10 cm above the column foot while ignoring the bond-slip effect between concrete and steel bars. Therefore, in order to improve calculation accuracy, two-scale models in this study are established using non-linear spring elements to consider the bond-slip effect in RC columns.

4. XFEM-Based Two-Scale Simulation on RC Columns under Cyclic Loading

4.1. Comparison of Failure Modes of Tested Specimens

The crack pattern of specimen C1 with ACR of 0.08 is shown in Figure 10a, presenting three transverse cracks at the heights of 5.0 cm, 15 cm and 30 cm, respectively, while loaded to 5 mm. With the increment of horizontal loading, oblique cracks appeared gradually and the highest height of the cracks was 45 cm from the bottom of specimen C1. The final failure pattern and crack distribution of specimen C2 with ACR of 0.10 are shown in Figure 10b. When loaded to 10 mm, oblique cracks occurred at the height of 15 cm from the bottom of specimen C2. When specimen C2 reached final failure mode, the highest height of the crack was 35 cm from the bottom of specimen C2. A total of three cross-cutting cracks appeared on the surface at loading sides, and the concrete at the column foot was crushed. The final failure pattern and crack distribution of specimen C3 with the ACR of 0.20 are shown in Figure 10c. Horizontal cracks occurred at the heights of 10 cm and 28 cm when loaded to 5.0 mm. When loaded to 20 mm, the concrete at the column root started crushing. As loading continued, the crack propagation was mainly concentrated at the bottom of specimen C3. It can be seen that the highest crack position in tested RC columns with greater ACR becomes lower and the concrete crushing will concentrate toward the column foot.

4.2. Time History of Strain and Crack Width

Taking the steel strain measured from a strain gauge located at 155 mm above the column foot of specimen C3 as an example, the comparison of time-history curves of strain and crack width obtained from numerical simulation and that observed from experimental measurement are presented in Figure 11; Figure 12. It can be seen that the strain obtained from the XFEM-based two-scale model is higher, especially for tensile strain. However, the peak values of the compressive strain are basically the same. As presented in Figure 12, the variation trends of the maximum crack widths in two-scale element model and solid element are consistent with each other. The relative errors between peak values of crack widths are less than 2%.

4.3. Relationship between Lateral Displacement and Reaction Force

Under lateral cyclic loading, the relationship between lateral reaction force and horizontal displacement at the loading point of tested specimens simulated with the solid element model, two-scale model and that of test results is summarized in Figure 13a–c. The peak values of reaction forces are summarized in

Table 3. Comparison of peak load values.

Model Type	Positive Peak Value (kN)	Positive Peak Error (%)	Negative Peak Value (kN)	Negative Peak Error (%)
C1-exprimental results	84.15	-	−87.66	-
C1-Two-scale	89.33	6.15	−93.87	7.09
C2-exprimental results	89.53	-	−93.06	-
C2-Two-scale	93.09	3.98	−97.33	4.58
C3-exprimental results	106.43	-	−102.33	-
C3-Two-scale	96.95	−8.90	−102.44	0.11

. It can be seen that the numerical simulation results are in good agreement with experimental results. The hysteresis loops corresponding to each loading cycle are consistent with that of tested specimens, presenting a remarkable pinch phenomenon. Figure 13d shows the comparison of the hysteretic behavior and the corresponding crack patterns of RC columns with different ACRs. From simulation results, four symmetrical cracks on both sides of the specimens C1 to C3 can be observed at the end of the cyclic loading. As shown in Figure 13d, the height of cracking position led by tension and compressive crushing area at the corner of column foot will be lower with the increment of the ACR, which is similar to experimental observation shown in Figure 10a–c. Therefore, the cracking pattern in tested RC columns can be accurately forecast using the XFEM-based two-scale modeling approach proposed in this study.

As listed in Table 3, the peak values of reaction forces in hysteresis curves obtained from the two-scale model are slightly higher than those of the solid element model. Certain errors exist between positive and negative peak values corresponding to the forward and reverse loadings. The fundamental reason for this phenomenon is that the fiber-beam element is incapable of simulating concrete cracking, and the reduction of material strength and stiffness can only be reflected according to the reduction of Young’s modulus and material strength in constitutive law. Since concrete is heterogeneous material, microcracking forms even at the initial loading stage due to high stress status at local regions, which cannot be simulated with fiber elements. Compared with the high-precision solid elements, a certain degree of error is inevitable. The maximum principal stress criterion is introduced to the material constitutive model of concrete to determine the initiation and propagation of cracking in XFEM elements. In ABAQUS, the asymptotic crack-tip function $F_{\alpha }\left(x\right)$ is neglected, which makes the crack-tip unable to stay inside the element and cut through at least one whole element once element stress reaches the maximum principal stress. Therefore, a slight difference in element stress may lead to the propagation paths slightly varying from each other.

Moreover, it should be noted that the peak value errors of reaction force between numerical simulation results and experimental data are within 9%, as illustrated in Table 3. The peak loads of a two-scale model of specimen C1 are 6.15% and 7.09% higher than that of experimental results. For specimen C3, the positive load peak of the two-scale model is 8.90% smaller than that of experimental results. The deviation between the FEM analysis and experimental observation is induced by multiple reasons, including the installation error of experimental equipment, initial defects in concrete, variation of boundary condition and incapability of simulating compressive failure of concrete of the XFEM method. More importantly, the random distribution of aggregates in concrete will make the tensile or compressive behavior and even the elastic Young’s modulus of concrete present certain discrepancies, as reported in the mesoscale numerical simulation performed by Chen et al. [3]. The influence of mesostructure variation of concrete on the macro-level structural behavior of RC columns will be discussed and compared in Section 5.

4.4. Computational Efficiency Analysis

The original intention of two-scale analysis is to improve computational efficiency with guaranteed calculation precision. In this study, all of the numerical studies were performed using Dell workstation equipped with 8 central processing units (CPUs) and 32 gigabytes memory. In hysteretic performance analysis simulated with solid element model and two-scale models, the element numbers and the corresponding CPU times of different models are compared in

Table 4. Comparison of element numbers and computation time of XFEM-based numerical models.

Specimen No.	Comparison Parameter	Solid Element Model (T1)	Two-Scale Model (T2)	Time Ratio (T1/T2)
C1	Number of elements	2258	1144	0.79
	Wall clock time(s)	3081	2444
C2	Number of elements	2258	1144	0.80
	Wall clock time(s)	3540	2848
C3	Number of elements	2258	1144	0.79
	Wall clock time(s)	3572	2825

.

The element number will be significantly reduced by using two-scale modeling method. However, the reduction of CPU time is not strictly proportional to that of element quantity since the main program of two-scale analysis has to call the external material subroutine developed in this study. The element quantities in a solid element model and two-scale model are 2258 and 1144, respectively. The ratios of CPU times between a two-scale model and solid element model corresponding to specimens C1 to C3 are 0.79, 0.8 and 0.79, which are closer to 0.80. Generally, the calculation speed can be improved by 20%. The time ratio corresponding to a CDP-based two-scale model and solid element model subjected to cyclic loading is 0.590, presenting higher computational efficiency. In solid element models and two-scale models of specimens C1–C3, the nodes of concrete elements and the longitudinal bars are connected using non-linear spring elements. It should be stated that the non-linear degree of numerical models and the element number will increase significantly while considering the bond-slip effect. Compared with simple RC columns investigated in this study, the improvement of calculation efficiency of the two-scale modeling approach will be more remarkable for large-scale RC frame structures.

5. Multiscale Simulation on RC Column under Monotonic Loading

5.1. Comparison of Failure Modes

The tested specimen subjected to monotonic loading was numbered as C4. The dimension and steel arrangement of specimen C4 were identical to that of C1, as detailed in Section 2. The ACR of C4 was also set to 0.08. The failure pattern of C4 is presented in Figure 14. Four cracks were clearly observed in the tested specimen C4, and the main crack was about 8 cm above the column foot. Severe concrete crushing occurred at the compressive side of column foot. In this section, the tested specimen C4 is constructed using the CDP model and XFEM method-based multiscale simulation strategies, respectively.

Figure 15 shows the tensile and compressive damage distribution patterns simulated using CDP-based multiscale models. It can be seen that the failure morphology at the tensile side and compressive crush basically agree well with that in experimental observations of the tested specimen C4. However, it is quite difficult to identify the main crack with the maximum width according to the tensile damage distribution presented in Figure 15. Moreover, it can be observed that the specific position of the element exhibiting tensile and compressive damage is different in multiscale models constructed with various mesoscale concrete samples, which denotes that local failure of RC columns is sensitive to the mesostructure variation of concrete.

Figure 16a–e show the comparison of crack patterns of different XFEM-based multiscale models subjected to monotonic loadings. For comparison, the solid element segments of XFEM-based two-scale models are constructed using homogeneous concrete material, and the corresponding crack patterns are presented in Figure 16f. Similar to the simulation results in Figure 15, the location and propagation direction of cracks simulated with XFEM-based multiscale models also change with the aggregate distribution of concrete at meso-level. Since the material strength of the ITZ is the weakest, the initial crack forms inside the ITZ and propagates along the outline of the aggregates.

Figure 16a–c are the crack patterns of multiscale models where numerical concrete has identical aggregate geometry but different distribution patterns, which indicates that the local position, length or even the extending direction of cracks are different due to the random distribution of aggregates. Meanwhile, the crack patterns in Figure 16d–e show that the geometric shape of aggregates also leads to a difference in the local crack patterns of RC columns. In contrast, the positions of prominent cracks in Figure 16c,d are closer to that of tested specimen C4 subjected to monotonic loading. The position of the main crack in Figure 16e corresponding to the multiscale model established with elliptical aggregates is obviously higher than that of practical experimental observations. The cracking morphology in a two-scale model shown in Figure 16f matches the best of that of tested specimen C4. However, the effect of random aggregate distribution on local behavior cannot be revealed using homogeneous material assumptions.

5.2. Relationship between Lateral Displacement and Reaction Force

In this section, the predicted load–displacement curves of tested specimen C4 simulated with multiscale models are discussed and compared with experimental observations, as summarized in Figure 17. The comparison of peak load values and the corresponding displacement is tabulated in

Table 5. Comparison of peak load values and corresponding displacement.

Model Type	Peak Value (kN)	Relative Error (%)	Displacement (mm)	Relative Error (%)
Tested specimen C4	99.7	-	19.82	-
CDP-based multiscale models
Circular sample 1	95.3	−4.4	17.25	−13.0
Circular sample 2	95.9	−3.8	16.42	−17.2
Circular sample 3	97.8	−1.9	18.04	−9.0
Elliptical sample	97.2	−2.5	18.52	−6.6
Polygonal sample	98.9	−0.8	21.22	−7.1
XFEM-based multiscale models
Homogeneous model	95.3	−4.4	16.29	−17.5
Elliptical sample	97.0	−2.7	23.28	17.5
Circular sample	99.9	0.2	25.22	27.2

. Figure 17a shows the load–displacement curves simulated with CDP-based multiscale models, presenting acceptable agreement with experimental results, especially in the elastic and ascending segments. In contrast, the peak values of reaction forces and the descending segments are more sensitive to the mesostructure variation of concrete. The maximum and minimum relative values of peak loads are −4.4% and −0.8%, respectively. The relative error of displacement corresponding to the peak loads between multiscale simulation results and test datum ranges from −6.6% to −17.2%.

Comparatively, the CDP-based multiscale models with polygonal aggregate samples match the best of the experimental results. Therefore, three sets of polygonal aggregate sample were established for XFEM-based multiscale models, and the corresponding load–displacement curves are illustrated in Figure 17b. It can be seen that the simulation results meet the test results at the ascending branch well. Moreover, the structural stiffness can be accurately predicted. However, the simulation encounters severe convergence problems. The convergence problems can be alleviated while using elliptical and circular aggregates in numerical concrete models, as shown in Figure 17c. In general, the load–displacement curves of XFEM-based multiscale models are very close to each other and to that of test results. The cracking point can be accurately predicted using the XFEM-based multiscale models constructed with elliptical and circular aggregates. As listed in Table 5, the relative errors of the peak reaction force between experimental observation and XFEM model constructed with homogeneous concrete, elliptical and circular aggregates are −4.4%, −2.7% and 0.2%, respectively. The relative errors in the displacement corresponding to load-bearing capacity in each model are −17.5%, 17.5% and 27.2%. Compared with CDP-based multiscale models, the XFEM-based multiscale models are more efficient in evaluating the global structural response and load-carrying capacity of RC columns.

Based on the multiscale analysis discussed above, it can be concluded that the cracking pattern and the load-bearing capacity of RC columns can be forecast with an acceptable accuracy using the XFEM-based multiscale modeling strategies developed in this study. The stochastic structural behavior of RC structures can be revealed from the inherent randomness in the mesostructure of concrete.

It has to be noted that 3D simulation of the whole RC specimen at meso-level is preferred since out-of-plane effects cannot be explicitly simulated in 2D models. However, the 3D mesoscale simulation is extraordinarily time-consuming and needs high configuration of a work station, which makes it impossible to perform with regular computers. Therefore, the 2D mesoscale simulation has been extensively employed in mesoscale simulation on RC structures. Besides, the feasibility and applicability of 2D mesoscale modeling have also been discussed in-depth [31]. For simplification, it is suggested to simulate RC column using plane stress elements [1]. According to the simulation results summarized in Section 5.1 and Section 5.2, it can be clearly seen that the 2D multiscale models are capable of predicting the failure pattern and structural response of tested RC columns. Compared with two-scale modeling method [16], the multiscale modeling method can predict the cracking evolution process accurately. Moreover, it is convenient to reflect on the fluctuation of structural responses of RC columns by using mesoscale concrete models with various aggregate samples

6. Concluding Remarks

In this study, the feasibility of an XFEM-based multiscale modeling approach was validated using experimental observations of four RC columns subjected to monotonic and hysteretic loadings. The influence of mesostructure variation in concrete on the macroscale behavior of tested RC column specimens was discussed in depth. The main conclusions of this study are as follows:

(1)

The simulation results show that ignorance of bond-slip effect in FEM models leads to relatively higher load-bearing capacity and structural stiffness of RC columns. Moreover, the specific position of the main crack in numerical models of RC columns is sensitive to bond-slip effect between concrete and steel bars.

(2)

The distribution pattern of main cracks in RC columns subjected to cyclic loading can be effectively predicted by XFEM-based two-scale models, which agrees well with the experimental results. Compared with solid element models, the computational efficiency of the XFEM-based two-scale model can be obviously improved with guaranteed simulation accuracy.

(3)

The influence of ACR on the structural behavior of tested RC column specimens is forecast utilizing the XFEM-based two-scale modeling approach. The numerical results match well with that observed from the experimental study. RC columns subjected to higher ACRs present cracks at lower parts of RC columns, greater structural stiffness and higher load-bearing capacity.

(4)

The monotonic behavior of RC columns can be simulated using XFEM-based multiscale models with acceptable accuracy. Compared with CDP-based multiscale models, the XFEM-based multiscale models are more efficient in predicting the cracking phenomenon and global structural response of RC columns. However, severer convergence problems also exist in XFEM-based multiscale simulations.

The XFEM-based multiscale modeling method provides an efficient way for the numerical simulation of the non-linear behavior of RC members due to the advantages of XFEM in the analysis of crack initiation and propagation in concrete. The feasibility of the XFEM-based multiscale modeling method in 3D numerical simulation and the difference between 2D and 3D models will be further investigated in the following studies. Moreover, the experimental measurement of material property for aggregates, the ITZ and mortar utilizing nanoindentation equipment will be the research focus of future work.