Investigation of over-nonlocal damage and interface cohesive models for simulating structural behaviors of composite UHPC-NC members

Investigation of over-nonlocal damage and interface cohesive models for simulating structural behaviors of composite UHPC-NC members

Structures 28 (2020) 2617–2632 Contents lists available at ScienceDirect Structures journal homepage: Investigat...

9MB Sizes 0 Downloads 0 Views

Structures 28 (2020) 2617–2632

Contents lists available at ScienceDirect

Structures journal homepage:

Investigation of over-nonlocal damage and interface cohesive models for simulating structural behaviors of composite UHPC-NC members Siqi Yuan a, b, Zhao Liu a, c, *, Teng Tong a, c a

School of Civil Engineering, Southeast University, Nanjing, China Department of Civil and Environmental Engineering, University of Maryland, College Park, MD, United States c Key Laboratory of Concrete and Prestressed Concrete Structures of Ministry of Education, Southeast University, Nanjing, China b

A R T I C L E I N F O Original content: com/s/1Fuctc9JJwnAmWGEtTke1OA Keywords: UHPC-NC Gradient-enhanced Over-nonlocal Cohesive


In this paper, a finite element (FE) modeling method is proposed to predict the structural behaviors of composite members of ultrahigh-performance concrete (UHPC) and normal concrete (NC). This FE modeling method in­ tegrates: (i) over-nonlocal gradient-enhanced damage model for UHPC and NC, capable of circumventing strong mesh dependency of the classical continuum damage model and spurious damage growth of the standard gradient-enhanced damage model; and (ii) the irreversible viscosity-regularized cohesive zone formulation to describe UHPC-NC interfacial behaviors. Linearization of relative governing equations is presented and the framework is realized on general FE platform ABAQUS through user subroutines UEL and UMAT. Three UHPCNC composite specimens were prepared to validate the versatility of the numerical method. Improved slant shear test was accompanied to realistically quantify the UHPC-NC bond strength. The height of UHPC substrate gov­ erned the overall behaviors and leaded to complete different failure patterns. It is validated that incorporating the material properties and bond behaviors in this paper predicts the composite specimens’ behaviors in an accurate and efficient way.

1. Introduction Ultra-high-performance concrete (UHPC) represents the ultra-high strength cementitious family exhibiting compressive strength above 150 MPa, tensile strength above 5 MPa, excellent post-cracking behav­ iors, low permeability and agreed durability, etc. [20]. Over the past decades, the superior mechanical and material properties render UHPC a promising material in modern infrastructures, in terms of high-rise buildings [47], highway and pedestrian bridges [12,22], pavements [54], tunnels [28], and nuclear reactor containments [29]. However, the expensive initial cost, which is around 3~4 times that of normal concrete (NC), hinders a wider application of full-depth UHPC members (Saleem et al., 2011). The UHPC-NC composite members provide a more cost-effective solution to utilize the excellent material properties of UHPC [63]. Such as NC-UHPC composite bridge deck with a UHPC substrate and a NC overlay, this composite system can take advantages of both NC’s compressive strength and UHPC’s agreed post peak behaviors for one hand, and for the other hand the cracking and maintenance cost could be apparently alleviated [64]. Several experi­ ments were performed on the structural behaviors of UHPC-NC com­ posite members, including the strengthening of “as-built” NC members

with UHPC overlay [41], and “newly-built” UHPC-NC members inten­ ded to accelerate the construction speed [23]. These composite members have been shown to achieve agreed structural behaviors. Of no doubt, UHPC-NC interface is a vital influencing factor to these composite members. Debonding or the fracture of UHPC in the interfa­ cial zone would lead to their premature failures [36]. Although UHPC is reported to work compatible with NC [49], disperse experimental ob­ servations were witnessed and a consistent conclusion cannot be drawn. Habel et al. [21] found that only a few failures occurred at the interface due to a stronger bond strength than the RC tensile strength. Similarly, only some horizontal microcracks without any interface debonding were observed among the experiments conducted by Al-Osta et al. [3]. Yin at al. [62] recorded that the beams strengthened with UHPC patch pre­ sented typical flexural failure and not any interface debonding was observed. In contrast, interface debonding was a common phenomenon among the specimens with UHPC serving as overlay [52,37]. On the other hand, finite element (FE) investigations for the com­ posite members still remains challenging [63]. By assuming a perfect bond at the interface, Sadouki et al. [43] simulated the UHPC-NC composite beam tests performed by Noshiravani and Bruhwiler [36] in the DIANA FEA software. Al-Osta et al. [3] and Safdar et al [44] also

* Corresponding author at: School of Civil Engineering, Southeast University, Nanjing, China. Received 19 May 2020; Received in revised form 13 October 2020; Accepted 15 October 2020 Available online 9 November 2020 2352-0124/© 2020 Published by Elsevier Ltd on behalf of Institution of Structural Engineers.

S. Yuan et al.

Structures 28 (2020) 2617–2632

assumed a perfect interface bonding behavior in their simulation work, which may lead to an overestimate on the ultimate strength of the composite members. Lampropoulos et al. [27] assumed a friction coef­ ficient of 1.5 and a cohesion of 1.9 MPa for the interface in ATENA software to model the concrete beams strengthened with a UHPC jacket. According to the author’s knowledge, most relevant numerical in­ vestigations rely on mature constitutive models embedded in various FE softwares. A main drawback is that limited model parameters are authorized to be input, and therewith these models, intended for NC, cannot fully represent UHPC’s mechanical behaviors, especially after cracking [65]. The same issue exists in the interfacial behaviors between UHPC and NC layers. On the other aspect, the classical continuum models are unable to provide meaningful post-peak results, e.g. the concrete damageplasticity model provided by the ABAQUS software [3] or the multilinear equations for concrete [44]. The strain-softening will cause loss of ellipticity which may lead to the unphysically meaningful solution on refinement of the spatial discretization [39,6]. To this end, numerous regularization techniques are proposed to obtain mesh-objective results, as cohesive zone model [42], crack band model [8], nonlocal-based damage model [39], phase-field model [25] etc. Cohesive formula­ tions require inserting of zero-thickness interfacial elements between different 2D or 3D solid elements. The micro-mechanisms of interfacial debonding and failure are therefore reflected as the cohesive traction, which is a function of local separation [16,5]. Despite the relatively easy implementation of the crack band model, this approach is only useful only if the fracture zone is deduced in advance. Furthermore, it is difficult to extent the model for mixed model failures [59]. Phase-field modelling is able to represent all stages of fracture, such as crack nucleation, initiation, propagation, and coalescence [53]. There emerges an explosion of scholarly interests in this method. However, there are still some obstacles for its mature applications in cementitious materials [35,24]. Nonlocal regularizing technique realizes non-locality in the constitutive model, either in integral forms [7,9] or in gradientenhanced damage models [48,11]. The implicit gradient form is numerically attractive since the governing equation can be discretized easily, and has demonstrated its ability to serve as the localization limiters during softening. Luzio and Baˇzant [30] pointed out the stan­ dard the standard nonlocal enhancement may not fully regularize the problem for certain material models. They adopted the weighted sun of the local and nonlocal values in the micro-plane model for concrete to achieve full regularization and circumvent the spurious damage growth [55] at high deformation level in the standard gradient-enhanced damage. The weight of nonlocal value is set greater than unity and therewith the term “over-nonlocal” [56]. Other approaches attaining the similar effects include strain-based transient approach and its simplified version [19,46], the damage model with evolving nonlocal interaction [33], a localizing gradient damage-based interaction [40], and the smoothing gradient damage model [34]. The main contribution of the present study is twofold. Numerically, a FE modeling method for UHPC-NC composite members is proposed, with fracture/softening behaviors of UHPC and NC, and their interface being taken into account in a robust and efficient way. This modeling method is implemented with the ABAQUS software through usersubroutines UEL and UMAT. Experimentally, three UHPC-NC compos­ ite specimens were tested and the effects of UHPC-NC bond strength were clearly witnessed. A novel slant shear test was accompanied and the realistic bond strength was obtained. This paper is organized as follows. Section 2 presented the funda­ mental issues for over-nonlocal gradient-enhanced damage model for UHPC and viscosity-regularized cohesive formulation for UHPC-NC interface. The linearization of respective governing equations and the FE implementation in ABAQUS through user subroutines UEL and UMAT are given in Section 3. A model parameter identification is also performed in this section. Section 4 describes the composite specimens and experimental setup. Finally, this proposed FE modeling method is

validated with respect to the experimental data in Section 5. It is note­ worthy that the FE framework is realized with the implicit algorithm, rendering it suitable for static or quasi-static problems. In addition, the source UEL code of the over-nonlocal gradient-enhanced damage model for the numerical illustration is provided for the educational purpose. 2. Constitutive models In order to simulate the simulate the structural behaviors of the composite UHPC-NC composite members, this paper selected overnonlocal gradient-enhanced damage model for UHPC and NC material. The viscosity-regularized cohesive model was used for simulating the interface between these two. The first step for this modeling is to determine the constitutive models for each part in UHPC-NC composite members. 2.1. Over-nonlocal Gradient-enhanced damage model for NC and UHPC The softening behavior of UHPC or NC is governed by a local scalar equivalent strain εeq . To solve the mesh-sensitivity problem, an intuitive approach is to replace the local strain εeq with a nonlocal one εeq over the whole domain Ω: ∫ 1 εeq = ∫ ψ (y, x) εeq dΩ (1) ψ (y, x) Ω Ω where ψ (y, x) is a weight function which is often assumed to be uniform. The gradient nonlocal formulation derives the governing equation for the nonlocal strain εeq of the strong form as [38]:

εeq − c⋅∇2 εeq = εeq in Ω


where ∇2 denotes the Laplacian operator; c denotes a positive gradient parameter of the dimension length squared. The following Neumann boundary condition is applied to guarantee that the average of εeq equals that of εeq , as: ∇εeq ⋅n = 0on ∂Ωt


Relying on the concept of equivalent strain [15], the Cauchy stress tensor σ is obtained through the following simple relationship:

σ = (1 − d)σ = (1 − d)E : ε


where σ is the effective stress tensor; d is a scalar damage variable (0⩽d⩽1); E is the fourth-order elastic stiffness tensor; and ε is the strain tensor. For strain-based damage models, the damage scalar d(κ) is a function of a history variable κ. The damage loading function f(⋅) is introduced and its relationship to the rate of history variable Δκ has to meet the Kuhn-Tucker loading-unloading condition: f (⋅)⩽0 Δκ⩾0 Δκf (⋅) = 0


The function f(⋅) is expressed as: f (̂ε eq , κ) = ̂ ε eq − κ


where the over-nonlocal scalar equivalent strain ̂ε eq is the weighted sum of nonlocal and local scalar strain, εeq and εeq , as: ̂ ε eq = m⋅εeq + (1 − m)⋅εeq


where the weight parameter m is chosen greater than unity to realize full regularization, therewith the term “over-nonlocal” [30,26]. Damage irreversibility is satisfied with the implementation of the damage loading function f(̂ε eq , κ) . The exponential form of damage variable d(κ) is adopted [31], which exhibit a long tail of post-peak stress–strain response, as: 2618

S. Yuan et al.

⎧ ⎨ d(κ) =

Structures 28 (2020) 2617–2632


κ ⎩ 1 − 0 [(1 − α) + αexp(− β(κ − κ0 ))] κ

κ⩽κ0 κ > κ0


ϕn = exp(1)σmax,0 δn

ϕt =

√̅̅̅̅̅̅̅̅̅̅̅̅̅ exp(1) τmax,0 δt 2


where σ max,0 and τmax,0 indicates respectively the interface normal and tangential strengths. The coupling parameters q and r for normal and tangential directions in Eq. (10) are governed by:

where κ0 is a damage threshold level; α is the residual strength of ma­ terial which means a small portion of load can still be carried for the fully-damaged material; and β mainly controls the slope of the softening curve. The uniaxial tensile stress–strain curves and their damage evo­ lutions with different values of β are characterized in Fig. 1. The proper expression of local scalar equivalent strain εeq remains an open issue (Nguyen et al., 2018), and the frequent methods include Mazar’s model, smooth Rankine method, and modified von Mises method [57]. In this framework, the last one is adopted for UHPC and NC, which yields the flowing expression: √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ k− 1 1 (k − 1)2 2 2k I1 + εeq = I + J2 (9) 2k(1 − 2v) 2k (1 − 2v)2 1 (1 + v)2

q = ϕn /ϕt

r = Δ*n /δn


where Δ∗n is the value attained by the normal displacement jump after complete shear separation whenTn = 0. In this paper, we focus on model II fracture, which is debonding of UHPC-NC interface. Therewith, we assume that the tractions in normal and tangential directions evolve separately. Two parameters, (ϕn ,σmax,0 ) and (ϕt ,τmax,0 ), determine the shape of the traction-separation curves in the normal and tangential directions, respectively. With potential ϕ at hand, the tangential traction Tn and Tt are obtained by taking derivation of the potential ϕ to the local displacement jump Δn and Δt , as: ( ( ) ) ∂ϕ ϕn ⋅Δn Δn ∂ϕ ϕt ⋅Δt Δ2t Tn = = exp − = = 2 exp − T (13) t ∂Δn δn ∂Δt δ2n δ2t δ2t

where I1 is the first invariant of strain tensor; J2 is the second invariant of deviatoric strain tensor; and v is the Poisson’s ratio. The parameter k is defined as the ratio of the compressive strength fc to the tensile strength ft of the material k=fc /ft . It should be noted that more sophisticate pointwise material models for cementitious materials has already been proposed continuously [14]. However, the purpose of this contribution is to propose a robust FE framework, rather than concentrating on a specific material model. The pointwise accurate NC models can be easily incorporate into the FE framework in the later work, as well as the constitutive model for UHPC which has not been developed effectively.

Consequently, the cohesive responses in the normal and tangential directions are illustrated in Fig. 2. As a rule of thumb, complete crack is attained when Δn ⩾5δn or Δt ⩾5δt , since cohesive tractions never equal to zero [13]. Cohesive models suffer from the convergence problems for post-peak analyses, as the occurrence of the elastic snap-back instability. To this end, the technique of viscous regularization proposed by Gao and Bower [18] is included in the study. A small additional viscous dissipation is added to the traction-separation laws, as ( ) ( ) d Δn d Δt Tn,v = Tn + ξn Tt,v = Tt + ξt (14) dt δt dt δt

2.2. Viscosity-regularized cohesive mode for UHPC-NC interface A cohesive zone formulation is employed to describe debonding of UHPC-NC interface. The pivotal ingredient of a cohesive model is the traction-separation law. In this contribution, the work by Xu and Nee­

where ξn and ξt are viscosity like parameters governing the viscous en­ ergy dissipation under normal and tangential loadings, respectively. The viscosity parameters own invalid physical energy dissipation process, but to regularize instabilities triggered by the initiation of a crack on the weak plane.

dleman [60] is adopted. The traction vector T = {Tn Tt }T would de­ cays exponentially with the local displacement jump vector

Δ = {Δn Δt }T . In a more detailed way, a potential energy ϕ is defined as: [ ( ){[ ] ] ( )} Δn Δn 1 − q r − q Δn Δ2 ϕ = ϕn +ϕn exp − − q+ 1− r+ exp − 2t r − 1 δn δn δn r − 1 δt

2.3. Nonlocal parameters The damage parameters cannot be considered as the identical ma­ terial constants without calibrating the nonlocal parameters m and c [66]. The gradient parameter c governs the range of nonlocal interaction and is related to the square of the constant internal length scale lc . The over-nonlocal weight parameter m could be considered as a numerical


where δn and δt are the characteristic opening lengths in normal and tangential directions, respectively, ϕn and ϕt denote the normal and tangential work of separation respectively, and are defined as,

Fig. 1. Schematic of the uniaxial tensile behaviors. 2619

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 2. Traction-separation law characterizing the cohesive zone.

value, and any value larger than unity could fully regularize the solution. To well understand the effect of parameters m and c, a bar of length 100 mm is taken for an illustration, see Fig. 3. The bar is subject to uniaxial tensile loading by prescribed displacement (0.08 mm) at two ends. The following damage parameters are assigned, withE0 = 40 GPa, v=0.2, κ0 =0.000075, β=300, α= 0.99, andk = 10. It is noteworthy that the Young’s modulus is deliberatively reduced by 5% in the middle 10 mm wide zone, to initiate the localized damage. The bar is respectively meshed with 40, 80 and 160 elements. Moreover, three different values of c withc = 1 mm2, 5 mm2, and 15 mm2, and two values of m, withm = 1.0 and 1.2, are selected. The first set of analyses embraces the bar with 80 elements andm = 1.0, but different values of c. Their damage profiles are plotted in Fig. 4a. Clearly, c=1 mm2 approaches the conventional local damage model, as the damage concentrates within the 10 mm middle zone. The nonlocality can be clearly witnessed for the results ofc = 5 mm2 and 15 mm2, as their damage zones overflow the 10 mm weaken zone. The effect of different c on the overall force–displacement relationship is illustrated as the normalized stress–strain curves in Fig. 4b. A smaller or bigger c would lead to the brittle or the ductile FE result, both would deviate from the real mechanical response. Mesh sensitivity is investigated by the second set of analyses, which embraces the bar discretized with different element size, but with identical gradient parameterc = 5 mm2, and over-nonlocal weight parameterm = 1.0. Profiles of the damage indicator d is plotted in Fig. 5a, and the normalized stress–strain curves are plotted in Fig. 5b. It is clearly that mesh-objective solution is obtained for the mechanical responses of quasi-brittle materials by the over-nonlocal gradientenhanced damage model. Effects of over-nonlocal weight parameter m are investigated by the third set of analyses, of which the bar is discretized with 80 elements and the identical gradient parameterc = 15 mm2 is assumed. The overnonlocal damage model seems to enhance the post-peak strength compared to the standard nonlocal damage model, see Fig. 6a. Never­ theless, this small difference can be solved by adjusting the values of β in Eq. (8). The benefit of the over-nonlocal damage model can be witnessed in terms of damage and strain distributions, see Fig. 6b and c respec­ tively. The over-nonlocal damage model suppresses the spurious

damage growth of the standard gradient-enhanced damage model. The effect is more obvious if we look at the strain profile in Fig. 6c. 3. Finite element implementation 3.1. Implementation aspects in ABAQUS The software ABAQUS is adopted for its powerful built-in nonlinear solver and automatic time-stepping schemes. The constitutive models for UHPC, NC and their interface are implemented through UEL subroutine. As shape functions for serendipity and cohesive elements are defined by users, ABAQUS cannot extrapolate variables from Gauss points to element nodes. Therewith, it is not possible to trivially display the re­ sults in ABAQUS/Viewer. To facilitate the post-processing, an auxiliary dummy mesh is adopted, consisting of standard ABAQUS elements with the identical nodes and integration points. In this auxiliary mesh, stress components are deliberatively set as zero through UMAT, so as not to influence the final solution. The data in the UEL for each time increment is stored as built-in array SVARS. At the termination of this increment, the data in SVARS is transferred into the built-in array STATEV for each corresponding element and integration point of the dummy mesh. The data transfer is realized by the COMMON BLOCK [32] in the UEL and UMAT subroutines. 3.2. Implementation of over-nonlocal gradient-enhanced damage model For the over-nonlocal gradient enhanced damage model, the set of governing equations in the strong form comprises of the equilibrium equation and the implicit gradient Helmholtz equation, as: ∇⋅σ + b = 0 in Ω


εeq − c⋅∇2 εeq = εeq in Ω


with the following Neumann and Dirichlet boundary conditions:

σ ⋅n = t∗ ∇εeq ⋅n = 0 on ∂Ωt


u(x, t) = u* ON ∂Ωu


where b represents an external load vector; t∗ and u* are prescribed force and displacement at the boundary. Therewith, the displacements u(x, t) and so-called nonlocal scalar equivalent strain εeq (x, t) remain as the two unknown fields over the domain Ω. The strong form equations are discretized using a Galerkin FE approach. The weak form equations are obtained through weighting Eqs. (15) and (16) respectively with wu and wε , and applying the

Fig. 3. Schematic of 1-D bar in tension. 2620

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 4. Simulation results of different c.

Fig. 5. Simulation results of different element size.

divergence theorem and the boundary conditions in Eqs. (17) and (18), as: ∫ ∫ ∫ ∇wTu : σ dΩ = wTu : b dΩ + wTu : t* dΓ (19)

linearize these equations to fit the consistent incremental-iterative Newton-Raphson solution method. At the element level, the lineariza­ tion for nodal variables at iteration i with respect to the previous iter­ ation i − 1 is given by:

u i = u i − 1 + δu ̃ ̃ ̃







∫ Ω

wTε εeq dΩ +


∫ ∇wTε c∇εeq dΩ =


wTε εeq dΩ


Discretization of the weak form equations at element level using linear element shape functions and their derivatives for the displace­ ment and nonlocal strain is formulated as: u = Nu : u and εeq = Nε : ε eq ̃ ̃


ε = Bu : u and ∇εeq = Bε ε eq




̃eq, i

and the stress tensor at the Gaussian point is:

σ i = σ i− 1 + δσ


After some manipulations, we can obtain: ⎡ ⎤ ∂ d ∂ κ ∂ε eq ⎦ :δε − m ∂d ∂κ E : εi− 1 :δεeq E: εi− 1 : δσ = ⎣(1− di− 1 )E− (1− m) ∂κ ∂̂ε eq ∂ε ∂κ ∂̂ε eq (29) Subsisting Eqs. (27)–(29) into Eq. (23), and replacing the variables at Gaussian point with the nodal variables, yields:


∫ (NTε : Nε + c⋅BTε : Bε ) : ε eq dΩ = Ω ̃


NTε : ε eq dΩ ̃


Remembering Eq. (4) characterizing the stress–strain relationship ε eq , the and Eq. (7) standing for the over-nonlocal equivalent strain ̂ increment of stress tensor δσ at the Gaussian point is expressed as: ( ) ∂d ∂κ mδεeq + (1 − m)δεeq E : εi− 1 δσ = (1 − di− 1 )E : δε − (28) ∂κ ∂̂ε eq

where u and ε respectively are the vectors of nodal displacement and ̃ ̃eq nodal nonlocal equivalent strain; Nu and Nε are the shape functions; Bu and Bε are their partial derivatives with respect to the spatial co­ ordinates. With respect to Bubnov-Galerkin method, the weight func­ tions wu and wε can be discretized analogously, yielding the following discretized equilibrium equations: ∫ ∫ ∫ BTu : σ dΩ = NuT : b dΩ + NuT : t* dΓ (23) Ω

=ε + δ ε eq ̃eq, i− 1 ̃


Eqs. (23) and (24) are nonlinear and therewith, we must further 2621

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 6. Simulation results of different over-nonlocal coefficient m.

⎧ ⎨∫


⎬ ∂d ∂κ ∂εeq ⎦ : Bu dΩ : δu + : ⎣(1 − di− 1 )E − (1 − m) E:ε: ⎭ ̃ ∂κ ∂̂ε eq ∂ε

κ⩽κ0 ⎨ 0 ∂d = κ0 κ0 ⎩ 2 [1 − α + αexp(− β(κ − κ0 ))] + ∂κ αβexp(− β(κ − κ0 ) κ > κ0

⎩ Ω ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏟ Kuu i− 1 ⎧ ⎫ ∫ ⎨ ⎬ ∂d ∂κ T − m B : E : εi− 1 : Bε dΩ : δεeq = ⎭ ⎩ Ω ∂κ ∂̂ε eq u


u fext

u fint, i−


The derivative of local scalar equivalent strain to the tonsorial strain (∂εeq /∂ε) can be refer to the appendix of [45] and is not presented for succinctness. Finally, the linearization of the governing equations at iteration i with respect to the previous iteration i − 1 is outlined as follows. ⎫ [ ]⎧ [ u ] [ u ] ⎬ f int, i− 1 f ext Kuu Kui−ε 1 ⎨ du i− 1 ̃ − (34) = ε ε ⎩ dε ⎭ f ext f int, Kεi−u 1 Kεε i− 1 i− 1 ̃ eq


Linearization of Eq. (24) is performed in a similar approach and is expressed as: {∫ ( } {∫ } ) ∂εeq − NεT : (NεT : N ε + c⋅BTε : Bε ) dΩ : δε̃ eq = : Bu dΩ : δu + ∂ε Ω Ω ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏟ ̃ ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏟ Kuu i− 1

Kεi−u 1

{∫ 0 − ⏟⏞⏞⏟ fint ε


The derivative of the history variable κ to the over-nonlocal scalar strain ̂ε eq depends on the status of loading or unloading, and is expressed as: { ∂κ 1 loading = (33) 0 unloading ∂̂ε eq

Kui−ε 1

∫ ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏟ ∫ ∫ NuT : b dΩ + NuT : t* dΓ − BTu : σi− 1 dΩ Ω Γ Ω ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅⏟ ⏟̅̅̅̅̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅̅̅̅̅⏟


} ∫ (NεT : N ε + c⋅BTε : Bε ) : ε̃ eq, i − 1dΩ − NεT : ε̃ eq, i − 1dΩ Ω Ω ⏟̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏞⏞̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅⏟

uε εu εε where the sub-matrices Kuu i− 1 , Ki− 1 , Ki− 1 and Ki− 1 , and the residual force u u ε ε vector f ext , f int, i− 1 ,f ext and f int, i− 1 are defined in Eqs. (30) and (31).

fext, i− 1 ε

(31) For the exponential softening Eq. (8), the derivative of the damage parameter to the history variable reads:

The sub-matrix Kui−ε 1 does not equal to the transpose of the sub-matrix

Ki− 1 , rendering the non-symmetry of the stiffness matrix. The quadratic rate of convergence is observed when the non-symmetric tangential stiffness matrix is used in a full Newton-Raphson algorithm [39]. In this contribution, we implement the over-nonlocal gradientenhanced damage model for 2D problems, and it is not difficult to extend εu


S. Yuan et al.

Structures 28 (2020) 2617–2632

{ tloc =

tt,loc tn,loc

} (39)

= Cloc : Δuloc

where Δuloc is the local separation; Cloc is the tangential stiffness given by the traction-separation law. The separation in the local coordinate system Δuloc is related to the separation in the global one Δu, as (40)

Δuloc = LT : Δu with


L = [tt , tn ]

where unit vectors tt and tn are the direction cosines of the local co­ ordinates to the global ones. The Δu at the Gaussian point is interpolated through the nodal separation in the global system, as: { } Δux (ξ) Δu = = H(ξ) : Δu (42) Δuy (ξ) ̃ where H(ξ) is a 2 × 6 matrix as: [ ] N1 (ξ) 0 N2 (ξ) 0 N3 (ξ) 0 H(ξ) = 0 N1 (ξ) 0 N2 (ξ) 0 N3 (ξ)

Fig. 7. Degrees of freedom of the 2D elements.

to 3D problems. The 2D 8-node serendipity elements are adopted. To achieve quadratic rate of convergence, the quadratic shape functions and their derivatives Nu and their derivatives Bu for the displacement field u(x, t) are defined over the eight node, see Fig. 7. On the other aspect, the linear Nε and Bε for the nonlocal strain field εeq (x, t) are defined over the four corner nodes. The nodal displacement vector u and ̃ nodal nonlocal equivalent strain vector ε eq are given as: ̃ [ ]T u = u1x , u1y , u2x , u2y , ...u8x , u8y (35) ̃

ε eq = [ε1eq , ε2eq , ε3eq , ...ε4eq ]

in which the Ni (ξ) is a standard shape function can be found in any FE textbook, and ξ( − 1⩽ξ⩽1) denotes the local coordinate of the element. The nodal separation Δu in the global system is derived from the ̃ nodal displacement vector u, as: ̃ Δu = [ − I6×6 ̃


N1 0

0 N1

... N8 ... 0

] [ N1 0 Nε = N8 0

0 N1


N4 0

and their derivatives with respect to spatial coordinates ⎤ ⎡ ∂N1 ∂N8 0 ... 0 ⎡ ⎥ ⎢ ∂x ∂x ∂N1 ⎥ ⎢ ⎥ ⎢ 0 ... ⎢ ∂x ⎢ ∂N1 ∂N8 ⎥ ⎢ ⎥ ⎢ ... 0 Bu = ⎢ 0 B = ⎢ ε ∂y ∂y ⎥ ⎣ ∂N1 ⎢ ⎥ 0 ... ⎢ ⎥ ∂y ⎣ ∂N1 ∂N1 ∂N8 ∂N8 ⎦ ... ∂y ∂x ∂y ∂x

0 N4




Combining (44) and (45) yields:

] (37)


Δu = H(ξ) : Φ : u = B(ξ) : u ̃ ̃


∂N4 ∂x


I6×6 ] : u = Φ : u ̃ ̃

u = [u1x , u1y , u2x , u2y , ...u6x , u6y ] ̃

with Nu =


where I6×6 is a 6 × 6 unity matrix; u is given automatically with the ̃ software as:





where B(ξ) is a 2 × 12 matrix transferring the nodal displacements to the separations at the Gaussian point. Consequently, the force vector and stiffness matrix at the element level, in the global coordinate system required by UEL in Abaqus, is obtained as: ∫ ∫ ∫1 ( T ) f el = B : L : tloc ⋅detJ⋅dξ BT : t dA = w BT : t dL = w

⎤ 0 ⎥ ⎥ ⎥ ∂N4 ⎦ ∂y



− 1

(38) We adopt nine Gaussian integration points for the two unknown fields. In numerical implementation of the above model, the solution procedure follows the standard incremental-iterative Newton-Raphson method with displacement-control loading scheme.


and ∫


BT :

Kel = w


) L : Cloc : LT : B⋅detJ⋅dξ

− 1


where w is the thickness of the cohesive elements.

3.3. Implementation of viscosity-regularized cohesive model

4. Description of experiments

To fit the eight-node 2D element for the bulk, the 2D quadratic line cohesive element is formulated to simulate the interface behavior, see Fig. 7. The zero-thickness cohesive element connects the surfaces of adjacent 8-node elements for UHPC substrate and NC overlay. The two surfaces of the cohesive element initially lie together in the unstressed deformation state (zero-thickness). The traction vector for }T { the cohesion element tloc = tt,loc tn,loc in the local coordinates is obtained as:

4.1. UHPC-NC composite specimens In order to further investigate the structural behavior of this inno­ vative UHPC-NC composite members, this paper fabricated three UHPCNC composite beam specimens. The rectangular cross-section of these specimens was 250 mm × 150 mm (width × height) and the whole length was 1200 mm. As shown in Fig. 8, the height of UHPC substrate in 2623

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 8. Illustrations of the UHPC-NC composite specimens (unit: mm).

specimens UN50, UN60 and UN70 were50 mm, 60 mm, and 70 mm respectively. To reveal the role of interface bond strength, it was note­ worthy that all these three specimens were not reinforced with any rebars. At first, the UHPC substrate was fabricated and cured under natural environmental for 28 days, see Fig. 9. In order to achieve a good interface bond behavior, all the UHPC substrate was submerged into water before casting the NC overlay [1]. This research is aim to reflect the role of interface on the structural behavior of UHPC-NC composite members. It is notable that not any roughness treatment was imple­ mented on the UHPC surface in this research [17]. After a 28-day natural curing, a force-control four-point bending test was performed on these UHPC-NC specimens. Thee clear span between the two steel roller supports was 1000 mm and the shear-span ratio is 1.5. Three traditional LVDTs sere set to evaluate the reflection at the supports and midspan. Besides, a digital image correlation (DIC) system was also equipped to monitor the crack propagation and present the strain field of the specimens during the test [51,50].

cylindrical strength offc = 45.9 MPa. The theoretical values of NC Young’s modulus and tensile strength for can be obtained from ACI [14] ′ as: ft = 5.5 MPa andEc = 31.8 GPa. The utilized UHPC is a commercial product provided by Sobute New Materials Co., Ltd. Nanjing, China. Mixture of one cubic UHPC embraced 2095 kg/m3 UHPC power, 156 kg/m3 20 mm-length steel fiber taking, 22.1 kg/m3 superplasticizer and 182.4 kg/m3 water. Nine 100 mm × 100 mm × 100 mm UHPC cubes were compressed and the mean cubic strength was 129.7 MPa, which was equivalent to the mean cylindric strength of 108.9 MPa. The 55.5 GPa elastic modulus was obtained for UHPC with three 1000 mm × 1000 mm × 300 mm prisms. Two 500 mm length custom-fabricated dog-bone specimens were fabricated to determine the UHPC’s tensile strength. This tensile specimen owned middle rectangular cross-section of 100 mm × 50 mm. The full stress–strain curves were obtained, see Fig. 10, characterizing the mean tensile strength of 5.52 MPa. The material properties are summarized in Table 1 below. According to the material property test mentioned in this section, the parameters values to govern the NC’s tensile and compressive behaviors are set as: Ec =31.8 GPa, v=0.2, κ0 = 0.00017, β=4500, α= 0.99, andk = 8.4. The parameters values to govern the UHPC’s tensile and compres­ sive behaviors are set as: Ec =55.5 GPa, v=0.2, κ0 = 0.0001, β=4500, α= 0.99, andk = 19.6. ′

4.2. Material properties The material properties of the UHPC and NC used in this research were also tested. The NC 28-day compressive strength is tested on nine 150 mm × 150 mm × 150 mm cubic samples. The mean value of the obtained cubic strength was 54.6 MPa, which was equivalent to the 2624

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 9. Prefabrication of composite specimens and experimental setup.

angles were considered, namely 50◦ , 60◦ and 70◦ . This geometry design can guarantee the interface failure, rather than the local NC damage recorded in the other researches. In addition, three layers of CFRP strip was wrapped around the upper part to avoid the NC compressive damage observed in the collected database [58]. Same with the UHPCNC composite members in the previous section, no roughness treatment on the surface was implemented. These slant shear specimens were naturally cured for 28 days. After that, a displacement-control loading was applied along the axles, from a 1000kN electrohydraulic servo machine. The magnitudes of the loading force and the vertical displacement at loading point were precisely monitored by the experiment system. The measured interfacial strength τtest in slant shear test can be ob­ tained as:

Fig. 10. Experimental and fitted behavior for UHPC tensile strength.

4.3. Slant shear test A variety of test methods can be used to obtain the bond strength of two materials, including: pull-off test, push-out test, splitting tensile test, and slant shear test. Compared with the other test methods, the charac­ teristic of slant shear test is that this method is able to access the interface bond strength under combined shear and compressive forces[4]. An improved slant shear test with nine specimens was conducted in this paper. The geometry details of these specimens are illustrated in Fig. 11. All these slant shear specimens owned 400 mm in height with a larger cross-section (150 mm × 100 mm) for lower part and a smaller one (100 mm × 100 mm) for upper part. Three different horizonal Table 1 Measured properties for NC and UHPC (unit: MPa). NC UHPC



Yong’s modulus

54.6 129.7

– 5.52

– 55.5

Fig. 11. Illustration of slant shear test (unit: mm). 2625

S. Yuan et al.

Structures 28 (2020) 2617–2632

The material properties for UHPC and NC are identical to these described in the previous sections. The nonlocal parameters are chosen as c = 20 mm2 and m = 1.2 for all the analyses. The three possible FE simulations are considered: (i) no bond between UHPC and NC layers, which means that only vertical contact is allowed; (ii) the realistic UHPC-NC bond behaviors derived from the slant shear test; (iii) perfect bond between UHPC and NC layers, which means no interfacial debonding is allowed. For the second case, interface properties listed in Section 4.3 are adopted. Implicit analysis which is suitable for static or quasi-static loading conditions was used to numerically solve the nonlinear problem. Displacement-control loading scheme is applied and the strain rate ef­ fect is not considered in this contribution.

Table 2 Slant shear strength results. Angle

70◦ 60◦ 50◦

Sample no.

Tested Peak load (kN)

Shear Stress (MPa)

1 2 3 1 2 3 1 2 3

351.1 416.8 243.2 29.82 488.7 463.5 109.9 556.7 492.3

9.93 12.16 6.78 1.30 14.09 13.36 4.67 18.27 16.16

Fitted stress (MPa)


8.72 9.99 6.64 2.72 14.65 13.99 5.20 18.40 16.49

1.14 1.22 1.02 0.48 0.96 0.95 0.90 0.99 0.98

5. Simulation results and verification

τtest =

P⋅sinα A0 /cosα

5.1. Specimen UN50


As shown in Fig. 13, failure pattern for specimen UN50 is a typical flexural failure with a main vertical crack being initiated and propagated between the two loading points. Of special attention is that no UHPC-NC interface debonding was observed during the test. As no reinforcing rebars were embedded in the tensile layer, the composite beam failed at the peak load of 70.7 kN. The first simulation deals with the FE model with the perfect bond condition, of which no debonding is allowed for the UHPC-NC interface. The damage distribution of the composite FE model is plotted in Fig. 14a. Alongside, the external load versus the mid-span deflection curve of the FE model is compared with tested one in Fig. 14b. It is found that the stiffness matches with each other, nevertheless, the peak load from the FE simulation is 106.7 kN, much bigger than the measured peak load of 70.7 kN. The main reason for the over-estimated strength is that homogenous material properties are assigned to these user-defined el­ ements. Therefore, two damage zones grow symmetrically in the leftside and the right-side zones. The FE results deviates from the observed phenomenon that only one main crack formed and propagated. Obviously, the two damage zones require much more external work than the single damage zone, resulting in a higher FE peak strength. Cementitious materials are the heterogeneous materials and their material properties are randomly distributed [61]. The above FE model cannot reflect the above characters. Therewith, stochastic FE analysis should be taken into account to supplement the proposed FE framework. This part of work will be incorporated in the further investigation and

where P is the peak force, and A0 is the contact are perpendicular to the load (A0 =100 mm × 100 mm). The measured P and τtest are summarized in Table 2. AASHTO [2] suggests that adhesion and friction attribute to the interfacial strength, as:

τ = c+μ

P⋅cosα A0 /cosα

= c + μ⋅σ


where c is adhesive bond, and μ is the friction coefficient. Least-square fitting method yields the following coefficients: c= 1.95 MPa andμ = 0.9. The corresponding C.V for the fitting result was 29.56%, and the correlation coefficient was 0.981, illustrating agreed experimental and fitting results, see Fig. 12. As a rule of thumb, c is within 1.7 ~ 2.8 MPa and μ is within 1.0 ~ 1.4 [2]. 4.4. Numerical models As shown in Fig. 7, the UHPC and NC layers are modelled with the 8node serendipity elements and their interface is captured by the 6-node quadratic line cohesive elements. A mesh of 10 mm is adopted for these models, which is compatible with the magnitude of the gradient parameter c.

Fig. 13. Final failure pattern for specimens UN50.

Fig. 12. Experimental and fitted UHPC-NC bond strength.


S. Yuan et al.

Structures 28 (2020) 2617–2632

herein, a relatively simple method is adopted. An artificial defect is introduced in the second FE simulation, by lowering the tensile strength of the one selected user-defined element, see Fig. 14a, to 90% of the target value. Alongside, the realistic bond properties are adopted for the UHPC-NC interface. The damage distribution of the second composite FE model is plotted in Fig. 14a, which agrees well with the experimental phenomenon. One major flexural-dominated crack occurs and no debonding of the UHPC-NC interface is found. Moreover, the external load versus the mid-span deflection curve for the FE model with the artificial defect matches well with the experimental curve, in terms of stiffness and peak load, see Fig. 14b. The FE peak load is 74.7 kN, which is compatible with the measured strength of 70.7 kN. The third FE model represents the situation without bond between UHPC and NC layers. Also introduced is the artificial defect in UHPC layer, as in the second FE mode. The corresponding damage distribution is plotted in Fig. 14a, of which the deformed configuration is plotted to illustrate the relative debonding between UHPC and NC layers, as marked in Fig. 14a. The external load versus the mid-span deflection curve is also comparatively plotted in Fig. 14b. The stiffness of the composite specimen without bonding decrease a lot, as well as the peak load of 52.7 kN. Now we focus on the damage evolution of specimen UN50. The second FE model with the artificial defect is selected for comparison, which yields relatively accurate simulation results. Damage evolution is experimentally observed with the assistance of the DIC system. The prepeak initiation and propagation of these cracks could not be observed by naked eyes. The flexural crack was initiated in the UHPC substrate when subject to the external load of 63.8 kN. The increasing load to 67.0 kN witnessed the cracking in the NC layer, whereas the crack inside the UHPC substrate remained relatively stable. when increased to 68.2 kN, these two cracks intertwined and leaded to the final major flexural crack. The FE damage evolutions according to the external load of 63.8 kN, 67.0 kN and 68.2 kN is comparatively plotted in Fig. 15, which agrees well with the experiment. 5.2. Specimen UN60 In specimen UN60, three features were observed: (i) the flexural crack were observed in the UHPC layer, however, it did not penetrate the over-all cross-section and the UHPC layer remained unbroken; (ii) the flexural crack was initiated in the NC layer and leaded to the breakage of the NC layer; and (iii) the UHPC-NC interface failed abruptly for the region marked in the figure, and the NC overlay detached totally from the UHPC substrate. It should be noted that the interfacial failure was the reason responsible for the NC overlay’s breakage. As no

Fig. 14. FE and experimental comparisons for specimens UN50.

Fig. 15. Damage evolution comparisons for specimens UN50. 2627

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 16. Final failure pattern for specimens UN60.

reinforcing rebars were embedded in the tensile layer, this UHPC-NC composite specimen failed at the peak load of 76.9 kN (Fig. 16). Benefit from the first FE simulation for the specimen UN50, we do not intent to perform the FE simulation without the artificial defect. Stochastic FE simulation should be further incorporated in the frame­ work to trigger the damage. As we observed the interfacial failure in the test, therewith, the first FE model deals with the situation with the perfect bond. An artificial defect is introduced in the FE simulation, see Fig. 17a. The FE damage distribution and external load-mid-span deflection relationship for the specimen UN60 are compared with the experimental data in Fig. 17b. If perfect bond is assumed, typical flexural failure can be reached and the FE peak load is 78.4 kN, which is close to the experimental peak load of 76.9 kN. In addition, the FE stiffness is also close to the experimental stiffness. The identical material properties and the artificial defect of the first FE model are assigned to the second FE model, except the interfacial behavior which adopts the experimental properties. The damage dis­ tribution and the load versus mid-span deformation curve from the FE simulation are plotted in Fig. 17a. The FE peak load is 74.8 kN, which is compatible with the measured strength of 76.9 kN. Nevertheless, completely different failure mode is observed. As marked in Fig. 17b, interfacial failure occurs when close to the peak load, rather than the merely flexural failure in the first FE model. The interfacial failure is also reflected in the load-mid-span deformation cure, of which a much more brittle post-peak descending curve is witnessed. No interfacial bond was assumed for the third FE model, but with UHPC/NC properties and artificial defect. The final damage distribution is plotted in Fig. 17a. The deformed configuration is plotted to reveal the relative slip between UHPC substrate and NC overlay. The relationship between the mid-span deflection and the external load is comparatively plotted in Fig. 17b. Without interfacial bond, we can clearly observe the decreased stiffness, and the reduced peak strength of 48.8 kN. We now turn to the second FE simulation which yields the most realistic numerical result to exploit this specimen’s damage evolution. Similar to specimen UN50, the DIC system is also employed to detect the crack initiation and propagation, especially these occurred before the peak load. The occurrence of the first flexural crack was witnessed in the UHPC substrate when subject to the external load of 72.0 kN. Mean­ while, the NC layer remained integrality without any cracks. Benefited from the UHPC’s agreed post-cracking behavior, this crack in the UHPC layer remained relatively stable when the external loading reached 74.0 kN. This load level triggered the flexural crack in the NC layer. With the increase of external load to its peak value of 76.9 kN, the flexural crack in the NC layer penetrate the cross-section. The underlying reason is the bond failure of UHNC-NC interface. Accordingly, damage evolutions from the FE simulation are comparatively plotted with the DIC images in Fig. 18, subject to the external load of 72.0 kN, 74.0 kN and 76.9 kN,

Fig. 17. FE and experimental comparisons for specimens UN60.

respectively. It is agreed that the FE simulation generally captures the realistic damage evolutions. It is noteworthy that the interfacial failure can be captured by the simulation, as marked in Fig. 18, which cannot be experimentally captured by the DIC system. 5.3. Specimen UN70 With the height of UHPC substrate increased to 70 mm, the following failure features was observed: (i) the shear failure occurred inside the shear span of the NC layers, rather than the flexural failure; (ii) the UHPC substrate remained overall integrity and no visible cracks were observed; and (iii) the UHPC-NC interface failed in the shear span, and the NC overlay detached from the UHPC substrate. It is reasonable to postulate that the interfacial failure leaded to the shear failure of the NC layer. Specimen UN70 failed at 100.5 kN, as no tensile rebars were provided in the tensile layer (Fig. 19). Similar to the specimen UN60, the first FE model for specimen UN70 2628

S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 18. Failure patterns and damage evolutions of specimen U60.

Fig. 19. Final failure pattern for specimens UN70.

embraces the realistic UHPC and NC properties, and perfect bond is assumed. The numerical damage distribution and the load versus deflection curve are compared with the experimental data in Fig. 20. The perfect bond assumption leads to the flexural failure of the UHPCNC composite beam at the peak load of 109.1 kN, which is close to the experimental peak load of 100.5 kN. In addition, the FE stiffness is also close to the experimental stiffness. The second FE model adopts identical UHPC and NC material properties, but the realistic interfacial bond behavior. The damage dis­ tribution and the load versus mid-span deformation curve from the FE simulation are plotted in Fig. 20. The FE peak load is 104.7 kN, which is compatible with the measured strength of 100.5 kN. Different from the first FE simulation, completely different failure pattern is witnessed. As marked in Fig. 20, interfacial failure occurs when close to the peak load, rather than the merely flexural failure in the first FE model. It is the interfacial failure that leads to the initiation of shear crack of the NC layer. Similarly, the third FE model assumes no interfacial bond between UHPC and NC layers. The final damage distribution is plotted in Fig. 20a under the deformed configuration. The relationship between the midspan deflection and the external load is comparatively plotted in Fig. 20b. Without interfacial bond, we can clearly observe the decreased stiffness, and the reduced peak strength of 55.9 kN. Damage evolutions of specimen UN70 from the second FE simulation are compared with the experimental data obtained from the DIC system.

Fig. 20. FE and experimental comparisons for specimens UN70.


S. Yuan et al.

Structures 28 (2020) 2617–2632

Fig. 21. Failure patterns and damage evolutions of specimen U70.

specimen UN50, the hybrid of flexural damage and interfacial failure for specimen UN60, and the hybrid of interfacial failure and shear crack in the NC overlay for specimen UN70. The interfacial behavior played a dominate role for the UHPC-NC composite members with a relatively thick UHPC substrate. (5) The proposed FE modeling method can realistically reflect the damage patterns and load–displacement curves of the three UHPC-NC composite members without rebars. Damage evolu­ tions predicted from the FE simulation are compared and vali­ dated with the experimental observations assisted by the DIC system. Artificial defect is introduced in the FE model to trigger the damage. Alongside, perfect bond, realistic bond and no bond are numerically investigated to demonstrate the role of the UHPC-NC interfacial behaviors.

The occurrence of some flexural micro-cracks inside the UHPC substrate was captured at the external load of 84.1 kN. The enhanced postcracking behavior of the UHPC rendered these flexural cracks stable, even until the termination of the test. The shear crack inside the NC overlay was observed at the external loading of 91.5 kN. Simulta­ neously, the NC-UHPC bond zone seemed to be not affected and no interfacial cracks were witnessed. When increased to 100.0 kN, the interface debonded and leaded to the abrupt propagation of the NC layer’s shear crack. It was interesting to find that with the UHPC height increased to 70 mm, the hybrid failure pattern of shear crack and interfacial debonding occurred, which was completely different from specimens UN50 and UN60. Accordingly, damage evolutions from the FE simulation are comparatively plotted with the DIC images in Fig. 21, subject to the external load of 84.1 kN, 91.5 kN and 100.5 kN, respec­ tively. It is agreed that the FE simulation satisfactorily reflect the un­ derlying failure mechanism of specimen UN70.

The versatility of the proposed FE framework is verified through the numerical and experimental investigations. Although a good bond quality is reported between UHPC and NC layers, it would not be equivalent to the perfectly bonded interface. It is demonstrated that the UHPC-NC interfacial behaviors are vital and to some extent would govern the composite member’s overall responses. Further work will incorporate UHPC-NC bond strength with different content of interfacial rebars, the UHPC-NC composite members with a high ratio of longitu­ dinal rebars, and the stochastic analysis of material properties. Besides, we admit that more detailed discusses on shrinkage and strain rate need to be conducted in the further research, which are not included in this manuscript.

6. Conclusion remarks The FE framework is proposed to model the structural response of composite UHPC-NC members, with the assistance of user subroutines UEL and UMAT in the ABAQUS software. To verify the versatility of the proposed FE model, four UNPC-NC composite specimens were loaded to failure. Through the numerical and experimental investigations, the following conclusions could be drawn: (1) The over-nonlocal gradient-enhanced damage model is employed for UHPC and NC materials, which can effectively regularize the mesh-sensitivity obstacle in the classical continuum damage model and spurious damage growth of standard gradientenhanced damage model. (2) The UHPC-NC interface is effectively and efficiently simulated with viscosity-regularized cohesive element. Merely two param­ eters are required to defined the interfacial behaviors. (3) The novel slant shear test was performed to obtain the real UHPCNC bond strength. The enlarged shape of NC part and the FRP warping were a good method to circumvent the typical NC failure in the literature. Therewith, accurate and realistic UHPC-NC bond strength was obtained. (4) Four UHPC-NC composite specimen were test. Three of them were not reinforced with any rebars in longitudinal or transverse directions. The increasing height of UHPC substrate witnessed three different failure patterns, namely pure flexural failure for

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. Data availability statement The source code, including INPUT file and user-subroutine UEL code for numerical illustrations in Sect. 5.3 is available to download from with pass­ word: z4fn.


S. Yuan et al.

Structures 28 (2020) 2617–2632


[30] Di Luzio G, Baˇzant ZP. Spectral analysis of localization in nonlocal and overnonlocal materials with softening plasticity or damage. Int J Solids Struct 2005;42 (23):6071–100. [31] Meschke G, Lackner R, Mang HA. An anisotropic elastoplastic-damage model for plain concrete. Int J Numer Meth Eng 1998;42(4):703–27. [32] Msekh MA, Sargado JM, Jamshidian M, Areias PM, Rabczuk T. Abaqus implementation of phase-field model for brittle fracture. Comput Mater Sci 2015; 96:472–84. [33] Nguyen GD. A damage model with evolving nonlocal interactions. Int J Solids Struct 2011;48(10):1544–59. [34] Nguyen THA, Bui TQ, Hirose S. Smoothing gradient damage model with evolving anisotropic nonlocal interactions tailored to low-order finite elements. Comput Methods Appl Mech Eng 2018;328:498–541. [35] Nguyen TT, Zhu YJ, Bornert M, Chateau C. A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure. Eng Fract Mech 2015;139:18–39. [36] Noshiravani T, Brühwiler E. Experimental investigation on reinforced ultra-highperformance fiber-reinforced concrete composite beams subjected to combined bending and shear. ACI Struct J 2013;110:251–61. [37] Paschalis SA, Lampropoulos AP, Tsioulou O. Experimental and numerical study of the performance of ultra high performance fiber reinforced concrete for the flexural strengthening of full scale reinforced concrete members. Constr Build Mater 2018;186:351–66. [38] Peerlings RHJ, de Borst R, Brekelmans WAM, Geers MGD. Gradient-enhanced damage modelling of concrete fracture. Mech Cohes-Frict Mater 1998;3(4):323–42. [39] Peerlings RHJ, Borst RDe, Brekelmans WAM, Vree JHPDe. Gradient enhanced damage for quasi-brittle materials. Int J Numer Meth Eng 1996; 39(19): 33913403. [40] Poh LH, Sun G. Localizing gradient damage model with decreasing interactions: localizing gradient damage model with decreasing interactions. Int J Numer Meth Engng 2017;110(6):503–22. [41] Prem PR, Murthy AR. Acoustic emission and flexural behaviour of RC beams strengthened with UHPC overlay. Constr Build Mater 2016;123:481–92. [42] Roesler J, Paulino GH, Park K, Gaedicke C. Concrete fracture prediction using bilinear softening. Cem Concr Compos 2007;29(4):300–12. [43] Sadouki H, Denari´ e E, Brühwiler E. Validation of a FEA model of structural response of RC-cantilever beams strengthened with a (R-) UHPFRC layer. Constr Build Mater 2017;140:100–8. [44] Safdar M, Matsumoto T, Kakuma Ko. Flexural behavior of reinforced concrete beams repaired with ultra-high performance fiber reinforced concrete (UHPFRC). Compos Struct 2016;157:448–60. [45] Sarkar S, Singh IV, Mishra BK, Shedbale AS, Poh LH. A comparative study and ABAQUS implementation of conventional and localizing gradient enhanced damage models. Finite Elem Anal Des 2019;160:1–31. [46] Saroukhani S, Vafadari R, Simone A. A simplified implementation of a gradientenhanced damage model with transient length scale effects. Comput Mech 2013;51 (6):899–909. [47] Schmidt M, Fehling E. Ultra-high-performance concrete: research, development and application in Europe. ACI Special Publ 2005;228:51–78. [48] Schreyer HL, Chen Z. One-dimensional softening with localization. J Appl Mech 1986; 791-797. [49] Shafieifar M, Farzad M, Azizinamini A. Experimental and numerical study on mechanical properties of Ultra High Performance Concrete (UHPC). Constr Build Mater 2017;156:402–11. [50] Shao X, Dai X, Chen Z, He X. Real-time 3D digital image correlation method and its application in human pulse monitoring. Appl Opt 2016;55(4):696. 10.1364/AO.55.000696. [51] Su Y, Zhang Q, Fang Z, Wang Y, Liu Y, Wu S. Elimination of systematic error in digital image correlation caused by intensity interpolation by introducing position randomness to subset points. Opt Lasers Eng 2019;114:60–75. [52] Tanarslan HM, Alver N, Jahangiri R, Yalçınkaya Ç, Yazıcı H. Flexural strengthening of RC beams using UHPFRC laminates: bonding techniques and rebar addition. Constr Build Mater 2017;155:45–55. [53] Tann´ e E, Li T, Bourdin B, Marigo J-J, Maurini C. Crack nucleation in variational phase-field models of brittle fracture. J Mech Phys Solids 2018;110:80–99. [54] Valipour M, Khayat KH. Coupled effect of shrinkage-mitigating admixtures and saturated lightweight sand on shrinkage of UHPC for overlay applications. Constr Build Mater 2018;184:320–9. [55] Vandoren B, Simone A. Modeling and simulation of quasi-brittle failure with continuous anisotropic stress-based gradient-enhanced damage models. Comput Methods Appl Mech Eng 2018;332:644–85. [56] Vermeer PA, Brinkgreve RBJ. A new effective non-local strain-measure for softening plasticity. International workshop on localisation and bifurcation theory for soils and rocks 1994. [57] de Vree JHP, Brekelmans WAM, van Gils MAJ. Comparison of nonlocal approaches in continuum damage mechanics. Comput Struct 1995;55(4):581–8. [58] Wei Y, Zhang Xi, Wu G, Zhou Y. Behaviour of concrete confined by both steel spirals and fiber-reinforced polymer under axial load. Compos Struct 2018;192: 577–91. [59] Xie De, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Fract Mech 2006;73(13):1783–96. [60] Xu X-P, Needleman A. Void nucleation by inclusion debonding in a crystal matrix. Modelling Simul Mater Sci Eng 1993;1(2):111–32. [61] Ye G, van Breugel K, Fraaij ALA. Three-dimensional microstructure analysis of numerically simulated cementitious materials. Cem Concr Res 2003;33(2):215–22.

This research is supported by the National Natural Science Founda­ tion (51778137). It is also supported by Shanghai Engineering Research Center of High Performance Composite Bridges (19DZ2254200). The support provided by China Scholarship Council (CSC) during a visit of the first author to University of Maryland, College Park is also acknowledged. References [1] Aaleti S, Sritharan S. Quantifying bonding characteristics between UHPC and normal-strength concrete for bridge deck application. J Bridge Eng 2019;24(6): 04019041. [2] AASHTO LRFD Highway Bridge Design Specifications American Association of State Highway and Transportation Officials 2007 Washington D.C. [3] Al-Osta MA, Isa MN, Baluch MH, Rahman MK. Flexural behavior of reinforced concrete beams strengthened with ultra-high performance fiber reinforced concrete. Constr Build Mater 2017;134:279–96. [4] ASTM-C822. Standard test method for bond strength of epoxy-resin system used with concrete by slant shear. 1999; West Conshohocken. [5] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7(1):55–129. [6] Bazant ZP. Instability, ductility, and size effect in strain-softening concrete. ASCE J Eng Mech Div 1976;102(2):331–44. [7] Baˇzant ZP, Belytschko TB, Chang T-P. Continuum theory for strain-softening. J Eng Mech 1984;110(12):1666–92. [8] Baˇzant ZP, Oh BH. Crack band theory for fracture of concrete. Mat Constr 1983;16 (3):155–77. [9] Baˇzant ZP, Pijaudier-Cabot G. Nonlocal damage: continuum model and localization instability. J Appl Mech 1988. [10] Blais PY, Couture M. Precast, Prestressed Pedestrian Bridge World’s First Reactive Powder Concrete Structure. pcij 1999;44(5):60–71. [11] Buildings code requirements for structural concrete (ACI 318). American Concrete Institute (ACI), 1999. [12] De Borst R, Mühlhaus H-B. Gradient-dependent plasticity: formulation and algorithmic aspects. Int J Numer Meth Engng 1992;35(3):521–39. [13] del Busto S, Beteg´ on C, Martínez-Pa˜ neda E. A cohesive zone framework for environmentally assisted fatigue. Eng Fract Mech 2017;185:210–26. [14] Dolado JS, van Breugel K. Recent advances in modeling for cementitious materials. Cem Concr Res 2011;41(7):711–26. [15] Cordebois JP, Sidoroff F. Damage induced elastic anisotropy. Dordrecht: Mechanical Behavior of Anisotropic Solids/Comportment M´ echanique des Solides Anisotropes. Springer; 1982. p. 761–74. [16] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8(2): 100–4. [17] Fehling E, Schmidt M, Walraven J, Leutbecher T, Frohlich S. Ultra-high performance concrete UHPC: Fundamentals, design, examples. 2014; John Wiley & Sons. [18] Gao YF, Bower AF. A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces. Modelling Simul Mater Sci Eng 2004;12(3):453–63. [19] Geers MGD, de Borst R, Brekelmans WAM, Peerlings RHJ. Strain-based transientgradient damage model for failure analyses. Comput Methods Appl Mech Eng 1998;160(1-2):133–53. [20] Graybeal B. Ultra-high performance concrete. No. FHWA-HRT-11-038 2011. [21] Habel K, Denari´e E, Brühwiler E. Experimental investigation of composite ultrahigh-performance fiber-reinforced concrete and conventional concrete members. ACI Struct J 2007;104(1):93. [22] Hajar Z, Simon A, Lecointre D. Construction of the first road bridges made of UHPC. 3rd International Symposium on HPC, Orlando 2003. [23] Honarvar E, Sritharan S, Matthews Rouse J, Aaleti S. Bridge decks with precast UHPC waffle panels: a field evaluation and design optimization. J Bridge Eng 2016; 21(1):04015030. [24] Hou Y, Guo M, Ge Z, Sun W, Wang L. Mixed-mode I–II cracking characterization of mortar using phase-field method. J Eng Mech 2017;143(7):04017033. https://doi. org/10.1061/(ASCE)EM.1943-7889.0001228. [25] Idesman AV, Levitas VI, Preston DL, Cho J-Y. Finite element simulations of martensitic phase transitions and microstructures based on a strain softening model. J Mech Phys Solids 2005;53(3):495–523. [26] Jir´ asek M, Rolshoven S, Grassl P. Size effect on fracture energy induced by nonlocality. Int J Numer Anal Meth Geomech 2004;28(78):653–70. [27] Lampropoulos AP, Paschalis SA, Tsioulou OT, Dritsos SE. Strengthening of reinforced concrete beams using ultra high performance fibre reinforced concrete (UHPFRC). Eng Struct 2016;106:370–84. [28] Li Ye, Tan KH, Yang E-H. Influence of aggregate size and inclusion of polypropylene and steel fibers on the hot permeability of ultra-high performance concrete (UHPC) at elevated temperature. Constr Build Mater 2018;169:629–37. [29] Lin Y, Yan J, Wang Z, Fan F, Zou C. Effect of silica fumes on fluidity of UHPC: experiments, influence mechanism and evaluation methods. Constr Build Mater 2019;210:451–60.


S. Yuan et al.

Structures 28 (2020) 2617–2632

[62] Yin H, Teo W, Shirai K. Experimental investigation on the behaviour of reinforced concrete slabs strengthened with ultra-high performance concrete. Constr Build Mater 2017;155:463–74. [63] Yin H, Shirai K, Teo W. Numerical model for predicting the structural response of composite UHPC–concrete members considering the bond strength at the interface. Compos Struct 2019;215:185–97.

[64] Yuan S, Liu Z, Tong T, Fu CC. A pilot study on structural responses of normal concrete-UHPC composite bridge decks w/w0 rebars through an experimentalnumerical approach. Adv Civ Eng 2020. [65] Zhu Y, Zhang Y, Hussein HH, Chen G. Numerical modeling for damaged reinforced concrete slab strengthened by ultra-high performance concrete (UHPC) layer. Eng Struct 2020;209:110031. [66] Zreid I, Kaliske M. A gradient enhanced plasticity–damage microplane model for concrete. Comput Mech 2018;62(5):1239–57.