Mechanical integrity and sustainability of pre-stressed concrete bridge girders repaired by epoxy injection

GEORGIA DOT RESEARCH PROJECT 17-08 FINAL REPORT
MECHANICAL INTEGRITY AND SUSTAINABILITY OF PRE-STRESSED CONCRETE
BRIDGE GIRDERS REPAIRED BY EPOXY INJECTION PHASE II
OFFICE OF PERFORMANCE-BASED MANAGEMENT AND RESEARCH
600 WEST PEACHTREE ST NW, ATLANTA, GA 30308

TECHNICAL REPORT DOCUMENTATION PAGE

1. Report No.: FHWA-GA-21-1708

2. Government Accession No.:

3. Recipient's Catalog No.:

4. Title and Subtitle: Mechanical integrity and sustainability of pre-stressed concrete bridge girders repaired by epoxy injection - Phase II

5. Report Date: December, 2020
6. Performing Organization Code:

7. Author(s): C. Arson, Ph.D. and K. Ji, M.Sc.

8. Performing Organ. Report No.: 17-08

9. Performing Organization Name and Address: Georgia Institute of Technology School of Civil and Environmental Engineering 790 Atlantic Drive NW Atlanta, GA 30332

10. Work Unit No.:
11. Contract or Grant No.: 0015712

12. Sponsoring Agency Name and Address: Georgia Department of Transportation Office of Performance-based Management and Research 600 West Peachtree St NW Atlanta, GA 30308

13. Type of Report and Period Covered: Final; October 2017 December 2020
14. Sponsoring Agency Code:

15. Supplementary Notes: Prepared in cooperation with the U.S. Department of Transportation and Federal Highway Administration

16. Abstract: During pre-stress transfer, the bottom portion of steel-reinforced concrete girders is subjected to an important compression induced by the relaxation of tension in the bars, which opens longitudinal cracks along the horizontal axis. Reactions at the supports induce shear stress, which sometimes translates into additional diagonal cracks at the ends of the girders. During the subsequent lifespan of the girder, a variety of crack patterns can occur, including longitudinal (along the beam axis), transverse (perpendicular to the beam axis), and diagonal cracks. At present, there is a need to assess the mechanical integrity and sustainability of pre-stressed concrete beams during the entire life cycle of the built infrastructure, including crack propagation and reparation. As such, this three-phase research program seeks to develop modeling strategies to predict the behavior of cracked concrete repaired by epoxy. The subject of this technical report is the Phase II research effort, which includes the numerical simulation of fracture patterns in concrete and the numerical assessment of damaged concrete reparation by epoxy injection. Preliminary models were presented in Phase I and the experimental data supporting the calibration and validation of the models proposed in Phase II is presented in the final report of Phase III. This report summarizes simulation results obtained at molecular, aggregate and metric scales. We show that HMWM deposition allows mechanical recovery of plain concrete, reinforced concrete and pre-stressed concrete members. Design recommendations are provided to optimize the use of HMWM for concrete reparation.

17. Key Words: Concrete, reinforced concrete, pre-stressed concrete girder, mechanical damage, fracture pattern, epoxy repair, simulation, design recommendations

19. Security Class (this report): Unclassified
Form DOT 1700.7 (8-69)

20. Security Class (this page): Unclassified

18. Distribution Statement: No restrictions

21. Number of Pages: 155

22. Price: Free

ii

GDOT Research Project 17-08
Final Report
MECHANICAL INTEGRITY AND SUSTAINABILITY OF PRE-STRESSED CONCRETE BRIDGE GIRDERS REPAIRED BY EPOXY INJECTION PHASE II
By
Chloe Arson, PhD, Associate Professor Koochul Ji, MSc., Ph.D. candidate
Georgia Tech Research Corporation Atlanta, Georgia
Contract with Georgia Department of Transportation
In cooperation with U.S. Department of Transportation Federal Highway Administration
December, 2020
The contents of this report reflect the views of the authors who are responsible for the facts and the accuracy of the data presented herein. The contents do not necessarily reflect the official views or policies of the Georgia Department of Transportation or the Federal Highway Administration. This report does not constitute a standard, specification, or regulation.
iii

TABLE OF CONTENTS

Page

LIST OF TABLES

VI

LIST OF FIGURES

VIII

EXECUTIVE SUMMARY

1

ACKNOWLEDGMENTS

6

CHAPTER 1. INTRODUCTION AND BACKGROUND

7

1.1 Introduction

7

1.2 Prior Phases

8

1.3 Phase II Objectives and Summary

10

1.4 Report Organization

11

CHAPTER 2. NUMERICAL SIMULATION OF FRACTURE PROPAGATION IN

PLAIN CONCRETE

13

2.1 DEM Model: Role of coarse aggregates on concrete strength

13

2.1.1 Introduction to the DEM

14

2.1.2 Displacement softening contact bond model

15

2.1.3 Mortar DEM model

17

2.1.4 Concrete DEM model

24

2.1.5 Sensitivity analyses

31

2.1.5.1 Effect of aggregate shape

31

2.1.5.2 Effect of the Interfacial Transition Zone (ITZ) on concrete strength 33

2.1.5.3 Effect of aggregate tensile strength on concrete strength

36

2.1.6 Summary and conclusions

39

2.2 CDP Model: Fracture patterns in plain concrete

41

2.2.1 Introduction to the CDP model

41

2.3 Verification of the CDP model against plain concrete experiments

46

2.3.1 Uniaxial compression tests on plain concrete cylinders

46

2.3.2 Notched three-point bending test without reinforcement

48

2.3.3 Three-point bending test on reinforced concrete beams

50

2.3.4 Summary and conclusions

53

CHAPTER 3. MULTI-SCALE STUDY OF PLAIN CONCRETE REPAIRED BY

EPOXY DEPOSITION

54

3.1 MD model of concrete/epoxy interface

54

3.1.1 Introduction to MD

54

iv

3.1.2 MD model of an interface of silica and highly cross-linked epoxy resin

(DGEBA and DETA)

61

3.1.2.1 Simulation of cross-linking and equilibration

61

3.1.2.2 Calibration of the MD epoxy model

64

3.1.2.3 Simulation of the silica/epoxy pull out test

70

3.1.3 Molecular models for the components of the concrete/HMWM interface 73

3.1.3.1 Molecular models of HMWM and concrete minerals

75

3.1.3.2 Material properties of HMWM

77

3.1.4 Molecular model of the mineral-HMWM interface

83

3.1.4.1 Construction of the MD interface models

83

3.1.4.2 Results of the numerical interface debonding tests with different mineral

substrates

84

3.1.4.3 Results of the numerical interface debonding tests under different

temperature and moisture conditions

88

3.1.4.4 Results of the numerical interface debonding tests under different strain

rates

90

3.1.4.5 Interface strength and fracture energy

91

3.1.4.6 MD predictions at low strain rate with phenomenological models

93

3.1.4.7 Contribution of the non-bonded energy to the work of separation of the

HMWM/mineral interfaces

96

3.1.5 Summary and conclusions

97

3.2 Simulation of concrete/HMWM at metric scale

98

3.2.1 HMWM model for FE simulation

98

3.2.2 Verification of the CDP model for concrete repaired by epoxy

100

3.2.2.1 Uniaxial compression tests on repaired concrete cylinders

102

3.2.2.2 Notched repaired three-point bending tests without reinforcement 105

3.2.3 Summary and conclusions

107

CHAPTER 4. NUMERICAL SIMULATION OF REPAIRED RC AND PSC BEAMS 108

4.1

Three-point bending tests with RC beams

108

4.1.1 Validation against experiments

109

4.1.2 Sensitivity analyzes

114

4.2 Numerical simulation of repaired prestressed concrete girder

117

4.2.1 Damage propagation during pre-stressing and loading

117

4.2.2 Simulation of pre-damaged PSC girders

124

4.2.3 Simulation of repaired PSC girders with pre-stress damage

126

v

4.2.4 Simulation of repaired PSC girders with end zone damage

129

CHAPTER 5. SUMMARY AND CONCLUSIONS

130

CHAPTER 6. REFERENCES

137

vi

LIST OF TABLES

Table 1. Mortar DEM model parameters

20

Table 2. CDP input stress-strain curve in compression

43

Table 3. CDP input stress-strain curve in tension

45

Table 4. Concrete material properties in the CDP model

46

Table 5. Possible Ensembles in MD

59

Table 6. Dimensions of crosslinked epoxy. W, D, H: model width, depth and height. 62

Table 7. Dimensions of the specimen used to simulate the pullout test at the interface between silica and highly cross-linked epoxy. W, D, H: width, depth, height. 70

Table 8. Dimensions of the simulation cells. W, D, H: width, depth, height

83

Table 9. Material parameters used in the FEM model of repaired concrete

101

Table 10. Bond parameters of the strands

119

Table 11. Bond parameters of the reinforcing rods

120

vii

LIST OF FIGURES

Figure 1. Force-displacement relationships used in the displacement softening

contact bond model

16

Figure 2. Stress-strain curves of mortar in the Brazilian tests simulated by the DEM 20

Figure 3. Micro-cracks in mortar at different stages of the Brazilian tests

21

Figure 4. Stress-strain curves of mortar in the uniaxial compression tests simulated

by the DEM

22

Figure 5. Micro-cracks in mortar at different stages of the uniaxial compression tests 23

Figure 6. Generation of aggregate shapes in the DEM model

24

Figure 7. DEM concrete specimens used to simulate (a) Brazilian tests and (b)

Uniaxial compression tests

26

Figure 8. Stress-strain curves of concrete in the Brazilian tests simulated with the

DEM

28

Figure 9. Micro-cracks in concrete at different stages of the Brazilian tests

29

Figure 10. Stress-strain curves of concrete in the uniaxial compression tests

simulated with the DEM

29

Figure 11. Micro-cracks in concrete at different stages of the uniaxial compression

tests

30

Figure 12. Effect of aggregate shape on concrete strength in Brazilian tests

32

Figure 13. Effect of aggregate shape on concrete strength in uniaxial compression

tests

33

Figure 14. Stress-strain curves of concrete with different interface contact ratios in

the Brazilian tests

34

Figure 15. Stress-strain curves of concrete with different interface contact ratios in

the uniaxial compression tests

35

Figure 16. A quasi-linear relationship between BTS and contact ratio

35

Figure 17. A quasi-linear relationship between UCS and contact ratio

36

Figure 18. Stress-strain curves of concrete with different aggregate tensile strength

in the Brazilian tests

37

Figure 19. Stress-strain curves of concrete with different aggregate tensile strengths

during the uniaxial compression tests

38

viii

Figure 20. A quasi-linear relationship between BTS and coarse aggregate tensile

stregnth

38

Figure 21. A quasi-linear relationship between UCS and coarse aggregate tensile

strength

39

Figure 22. Principle of the CDP model built in ABAQUS (a) compressive plasticity

(b) tension plasticity

42

Figure 23. (a) Experimental setup of the uniaxial compression test (b) FEM model

to simulate the uniaxial compression tests.

46

Figure 24. Stress-strain curve obtained by simulating the uniaxial compression tests

with the CDP model (black dots) in comparison with the input constitutive law

(blue line) and with the experimental strength (red dashed line)

47

Figure 25. Damage distribution during the uniaxial compression test

48

Figure 26. Three-point bending tests on notched plain concrete beams. Top:

Experimental set-up. Bottom: numerical model

49

Figure 27. Distribution of tensile damage during the three-point bending test

simulations conducted with the CDP model under a vertical displacement of

0.288 mm (left) and 3.08 mm (right).

50

Figure 28. Numerical and experimental force-CMOD curves obtained during the

three-point bending tests. Left: Sensitivity of the model to the stress ratio at the

kink point. Right: Sensitivity of the model to the tensile strength

50

Figure 29. Geometry of the reinforced concrete beam. The rebars were placed at a distance of 45 mm from the bottom of the beam and the spacing was 60 mm. 51

Figure 30. Damage evolution during the three-point bending test conducted on

reinforced concrete beams. Top: Numerical results (under displacements of 1.9

mm and 4mm). Bottom: Photograph of the experiment showing a specimen at

failure

52

Figure 31. Force displacement curves obtained numerically and experimentally for

the three-point bending test conducted on reinforced concrete beams

52

Figure 32. Examples of interface interlocking (Sasse & Fiebrich, 1983)

55

Figure 33. Multiscale analysis of epoxy-silica interface: "rubber-like model" (Lau et

al., 2012)

56

Figure 34. MD simulation flow chart (Lau et al., 2012)

57

Figure 35. Epoxy before and after Crosslinking between EPON-8672 and DETDA

(Li and Strachan, 2010)

60

Figure 36. Snapshots of atomic configurations during pullout tests simulated with

MD for copper-epoxy interfaces (Shaorui Yang, 2010)

60

ix

Figure 37. Molecular structures of DGEBA (top) and DETA (bottom). Grey atoms

are carbon, white atoms are hydrogen, red atoms are oxygen, blue atoms are

nitrogen

61

Figure 38. Epoxy after crosslinking between DGEBA and DETA (Arab and Nejad,

2012)

63

Figure 39. MD specimen of highly cross-linked epoxy. D = W = 55 ; H = 110 64

Figure 40. Properties of 78% cross-linked epoxy during the annealing simulation.

(a) Density vs. Temperature (b) Volume change vs. Temperature

65

Figure 41. Simulation of a uniaxial tension test on 90% cross-linked epoxy at 300K,

atmospheric pressure and using a strain rate of 5 x 108 s-1. (a) Stress-strain

response (b) Lateral vs. longitudinal strain curve

67

Figure 42. Stress-strain response of the 90% cross-linked epoxy subjected to a

uniaxial tension test at 300K, under different strain rates

68

Figure 43. Stress-strain response of cross-linked epoxy under a uniaxial tension test

imposed at a rate of 5 x 108 s-1 under 300K, for different conversion ratios

69

Figure 44. Stress-strain response of cross-linked epoxy under a uniaxial tension test

imposed at a rate of 5 x 108 s-1 under 300K, for two different sample sizes

70

Figure 45. Geometry of the interface between highly cross-linked epoxy and silica 71

Figure 46. MD simulation results obtained after simulating a pullout test on an

interface between silica and 90% cross-linked epoxy. Top: Stress-strain curve.

Bottom: Atomic configurations at key points of the stress-strain curve.

72

Figure 47. Stress-strain curves obtained after simulating a pullout test on an interface

between silica and cross-linked epoxy, for different conversion ratios

73

Figure 48. Molecular structure of (a) Cumene hydroperoxide (CHP), (b) Methyl

methacrylate (resin)

74

Figure 49. Schematic diagram of the chemical reactions involved in HMWM

crosslinking

75

Figure 50. HMWM atomic configuration after crosslinking. Carbon atoms are grey,

oxygen atoms are red, hydrogen atoms are white

76

Figure 51. MD lattice models of representative concrete minerals.

77

Figure 52. Simulation of an annealing test on cross-linked HMWM, at atmospheric

pressure, decreasing the temperature from 600K to 300K. Top: Volumetric

deformation vs. Temperature. Bottom: Density vs. Temperature

79

Figure 53. Simulation of a uniaxial tension test on cross-linked HMWM at 300K, at

atmospheric pressure and using a strain rate of 5 x 108 s-1. Top: Stress-strain

response. Bottom: Lateral vs. longitudinal strain curve

81

x

Figure 54. Effect of strain rate on the mechanical behavior of HMWM.

82

Figure 55. Molecular interface systems for (a) the calcite-HMWM model (b) the

silica-HMWM model.

84

Figure 56. HMWM/mineral interface debonding tests simulated at 300K, at a rate

of 5 x 108 s-1, for two different concrete minerals silica and calcite.

85

Figure 57. Calcite-water-HMWM interface system during the interface debonding

tests at 300K, at a rate of 5 x 108 s-1, with 100 water molecules.

86

Figure 58. Silica-water-HMWM interface system during the interface debonding

tests at 300K, at a rate of 5 x 108 s-1, with 100 water molecules.

87

Figure 59. Interface strength predicted by MD under different temperatures.

88

Figure 60. Interface strength predicted by MD under different moisture conditions 89

Figure 61. Effect of strain rate on the mechanical behavior of the HMWM/mineral

interface during the interface debonding tests

90

Figure 62. Simulation result and curve fitted exponential traction-separation law (in

red) to the stress-strain response obtained during the interface debonding test

conducted with the dry HMWM/silica interface system.

92

Figure 63. Variation of the tensile modulus of HMWM as a function of strain rate

in a uniaxial tension test. Experimental data is used as reference.

94

Figure 64. Interfacial strength as a function of strain rate and temperature using

Johnson-Cook model for the HMWM/calcite interface model. The red dot is an

experimental data point. The blue dots are the data points obtained by MD

simulation

95

Figure 65. Energy distribution during the interface debonding tests at 300 K for a

strain rate of 5 x 108s-1 under the dry condition

96

Figure 66. HMWM linear softening traction-separation law

98

Figure 67. HMWM-repaired concrete specimens used in the uniaxial compression

tests

102

Figure 68. Comparison of numerical and experimental results for the UC test. (a)

Stress-displacement curves obtained numerically, and strengths obtained

experimentally. (b) Numerical strength of the repaired specimen for different

CZ thicknesses. (c) Numerical stress-displacement curves obtained with

repaired concrete specimens, for different values of the exponent . (d) Contour

plot of damage at the end of the UC test simulation in repaired concrete,

compared to a photograph of the fractures observed experimentally

104

Figure 69. Non-repaired (top) and repaired (bottom) non-reinforced concrete beams

after three-point bending tests

105

xi

Figure 70. Evolution of the Crack Mouth Opening Displacement during the TPB tests done on plain and repaired non-reinforced concrete beams: numerical predictions are representative of the average behavior observed experimentally. 107

Figure 71. Photographs of the RC beams during the TPB tests and the corresponding

field of maximum principal strain obtained by DIC (a) when major cracks

developed (loading up to 100 kN); (b) when the beam with repaired major

cracks failed; (c) when minor cracks developed (loading up to 80 kN); (d) when

the beam with repaired minor cracks failed

110

Figure 72. Method to identify the damaged zone after the TPB tests on RC beams.

Top: the field of maximum principal strain is calculated by DIC. Bottom: the

zone where the maximum principal strain exceeds 0.002 is binarized and

meshed with SolidWorks

111

Figure 73. Numerical and experimental load/displacement curves obtained during

TPB tests performed on (a) As-built RC beams; (b) Repaired RC beams with

major cracks; (c) Repaired RC beams with minor cracks

113

Figure 74. Experiment and simulation results obtained for repaired RC beams: (a)

Photograph of repaired RC beam with major cracks; (b) Photograph of repaired

RC beam with minor cracks; (c) Snapshot of the FEM mesh, showing the zones

of damage at failure for the repaired beam with major cracks; (d) Snapshot of

the FEM mesh, showing the zones of damage at failure for the repaired beam

with minor cracks; (e) DIC results at failure, showing the field of maximum

principal strain in the repaired beam with major cracks; (f) DIC results at

failure, showing the field of maximum principal strain in the repaired beam

with minor cracks; (g) Snapshot of the FEM mesh, showing the distribution of

maximum principal strain at failure in the repaired beam with major cracks; (h)

Snapshot of the FEM mesh, showing the distribution of maximum principal

strain at failure in the repaired beam with minor cracks

114

Figure 75. Simulation of the TPB test with RC beams, after repairing cracks of

different widths. (a) Maximum load 100kN, crack width threshold 0.01 mm;

(b) Maximum load 80kN, crack width threshold 0.01 mm; (c) Load-

displacement curve for RC beams loaded up to 100 kN, with repaired cracks of

width > 0.01 mm, 0.05 mm, 0.1mm; (d) Load-displacement curve for RC

beams loaded up to 80 kN, with repaired cracks of width > 0.005 mm, 0.01

mm.

116

Figure 76. Details of prestressed concrete test beam (mm) (O. Yapar et al. 2015) 118

Figure 77. Applied loading configuration of the test beam (mm)

118

Figure 78. Exponential traction separation law (mm) (O. Yapar et al. 2015)

119

Figure 79. Simulation result of prestressed concrete girder (a) stretched strand to initial prestressing force level (b) longitudinal endzone damages followed by releasing of the strands (c) stress behavior of concrete followed by releasing of

xii

the strands (d) stress behavior of rebar and strands followed by releasing of the

strands

121

Figure 80. Simulation of a concrete girder subjected to prestressing followed by

four-point bending. (a) Comparison of the force vs displacement curves

obtained numerically to that reported in Yapar's experiments. (b) Force vs

displacement curves with prestress and without prestress (FEM results). (c)

Transferred shear stress from strands to concrete. (d) Effective stress in the

strands.

134

Figure 81. End zone damage configurations described by the PCI (a) Diagonal crack

- dia (b) Flange bottom - bot (c) Vertical crack - ver (d) Longitudinal crack -

longi

124

Figure 82. Damage distribution after pre-stressing, strand release and four point

bending of PSC girders with initial damaged end zones. (a) Diagonal crack -

dia (b) Flange bottom crack - bot (c) Vertical crack - ver (d) Longitudinal crack

longi

125

Figure 83. Load-displacement curves of PSC girders without and with initial damage 125

Figure 84. Post-processing of simulation result on PSC girder (a) Damage after

releasing prestressed strands (b) Zone where crack width is at least 0.01 mm

after releasing prestressed strands (grey color) (c) Binary image using a

0.01mm crack width threshold (d) Concrete part with cracks less than 0.01 mm

in width

127

Figure 85. Load vs. displacement curve of as-built PSC girder and PSC girder with

initial pre-stressed damaged zone repaired by HMWM deposition

128

Figure 86. Load vs. displacement curve of as-built PSC girder, RC girder of same

rebar volume fraction and RC girder with initial pre-stressed damaged zone

repaired by HMWM deposition

128

Figure 87. Load-displacement curves of PSC girders without and with initial damage

repaired by HMWM deposition

129

Figure 88. Effect of aggregate strenght and interface transition zone contact ratio

on concrete ensile strenght and compressive strength

134

Figure 89. Effect of concrete mineral composition, temperature and moisture on

concrete/HMWM interface strength.

135

Figure 90. Effect of HMWM reparation on the load capacity of RC beams

135

Figure 91. Effect of HMWM reparation on the load capacity of PSC beams

135

xiii

EXECUTIVE SUMMARY
During pre-stress transfer, the bottom portion of steel-reinforced concrete girders is subjected to an important compression induced by the relaxation of tension in the bars, which opens longitudinal cracks along the horizontal axis. Reactions at the supports induce shear stress, which sometimes translates into additional diagonal cracks at the ends of the girders. During the subsequent lifespan of the girder, a variety of crack patterns can occur, including longitudinal (along the beam axis), transverse (perpendicular to the beam axis), and diagonal cracks. At present, there is a need to assess the mechanical integrity and sustainability of pre-stressed concrete beams during the entire life cycle of the built infrastructure, including crack propagation and reparation. As such, this three-phase research program seeks to develop modeling strategies to predict the behavior of cracked concrete repaired by epoxy.
The subject of this technical report is the Phase II research effort, which includes the numerical simulation of fracture patterns in concrete and the numerical assessment of damaged concrete reparation by epoxy injection. Preliminary models were presented in Phase I and the experimental data supporting the calibration and validation of the models proposed in Phase II is presented in Phase III final report.
First, mortar and concrete cylinders were modeled with the Discrete Element Method (DEM). Mortar DEM parameters were calibrated against Phase III mortar cylinder tests. Concrete aggregates were then scanned. Point clouds obtained by scanning were imported in MATLAB to reconstruct the shape of aggregates within the DEM mortar models, so as to match the aggregate size distribution of GDOT concrete. The DEM parameters of the aggregates and of the Interfacial Transition Zone (ITZ) between the aggregates and the
1

mortar were then calibrated against the results of concrete cylinder tests performed in Phase III. DEM analyses allowed us to understand the effects of aggregate mechanical properties and of ITZ adhesive properties on fracture propagation and failure mechanisms in the concrete used by the GDOT. Notably, it was found that aggregates used in the composition of the concrete act as weak inclusions that decrease concrete strength.
A Molecular Dynamics (MD) model was designed to simulate pull out tests conducted on a concrete/DEGBA-epoxy and concrete/High Molecular Weight Metacrylate (HMWM) interfaces. It was found that the mechanical properties of cross-linked epoxy highly depend on the conversion ratio (i.e. on the cross-linking rate), the size of the specimen and the loading rate. Simulations of uniaxial tension tests on concrete/HMWM specimens showed that all interfaces fail by debonding, at the surface of contact between HMWM and the mineral substrate, and that the interfacial strength decreases in the presence of moisture, under low strain rate, or at high temperature. Silica/HMWM interfaces were systematically stronger than the calcite/HMWM interfaces. The results showed that the work of separation was mostly done by van der Waals forces. We showed that our MD results could be extrapolated to predict the tensile modulus of HMWM at low strain rate. The Johnson-Cook model was used to predict the HMWM/mineral interfacial strength for a broad range of temperatures and strain rates. MD results confirm that HMWM should be applied on dry surfaces and in concrete exposed to mild temperatures. HMWM application in concrete with high silica contents is more likely to last than in concrete with high calcite contents.
The Concrete Damage Plasticity (CDP) model was used to simulate fracture propagation in concrete at the metric scale. The model was calibrated against experimental
2

results obtained in Phase III for plain concrete cylinder tests and plain unreinforced concrete three-point bending tests. We enriched the CDP Finite Element model with Cohesive Zone (CZ) elements to represent the adhesion strength at the faces of cracks filled with epoxy. The model was calibrated against uniaxial compression tests and three-point bending tests performed on cut and repaired specimens in Phase III.
We used Digital Image Correlation (DIC) to identify the zones of high maximum principal strain after the first loading cycle applied to RC beams tested in Phase III. These zones were assigned repaired concrete elements and HMWM CZ elements for simulating the second load cycle, after reparation. The distribution of damage and the field of maximum principal strain calculated numerically at failure match the damage zone and the maximum principal strain field found experimentally. Simulation results suggest that HMWM can penetrate cracks of width 0.01 mm and above by gravity. We also found that HMWM reparation increases concrete stiffness and strength if cracks in concrete members are over 0.1 mm in width, in which case, the load capacity of repaired RC beams is 30 to 40% higher than that of as-built RC beams.
Models of Pre-stressed Steel reinforced Concrete (PSC) girders were built with the Finite Element Method (FEM) to simulate pre-stressing, strand release and four-point loading. The FEM model matches results reported in the literature and captures the development of damage during pre-stressing. We found that the load capacity of a PSC girder damaged by pre-stressing and then repaired would be about 7% higher than that of the as-built PSC girder. At same volume fraction of rebar, a Reinforced Concrete (RC) girder has a load capacity that is about half of that of the PSC girder. Addition of HMWM in a zone of damage similar to the pre-stressing damage zone would increase the load
3

capacity of the RC girder by 50%, reaching 75% of the load capacity of the as-built PSC girder. We simulated pre-stressing followed by strand release and four-point bending in PSC girders with initial damaged end zones, as described in the PCI guidelines. Compared to the as-built girders, pre-damaged girders exhibited a loss of stiffness of 20 to 25 % and a loss of load capacity of up to 15%. Simulation of repair followed by reloading shows that full recovery of stiffness and load capacity by HMWM deposition is possible.
Based on the results above, we suggest the following design recommendations: 1. Improve concrete
Use higher aggregate strength Use smaller aggregates Increase silica contents 2. Optimize the HMWM deposition technique Apply at low temperature Apply on dry crack faces Repair cracks > 0.01 mm in width 3. Do repair PSC beams according to PCI guidance Mechanical properties of PSC beams can be recovered upon reparation of
cracks classified by the PCI Reparation has a greater influence on load-displacement curve for lower
rebar volume fraction and lower prestress Load capacity of repaired RC beams with no prestress can exceed that of
as-built RC beams
4

4. Plot strength charts for in-situ implementation Full scale testing is needed to validate repaired PSC beam simulations Full scale parametric studies are needed to assess a larger set of crack distributions and reinforcement/pre-stress designs
5

ACKNOWLEDGEMENTS The following individuals at the GDOT provided many valuable suggestions throughout this study: Mr. David Jared, Mr. Christopher Watson and Dr. Peter Wu. The opinions and conclusions expressed herein are those of the authors and do not represent the opinions, conclusions, policies, standards, or specifications of the GDOT or of other cooperating organizations. The numerical research described herein was performed at the Georgia Institute of Technology Damage Poro-Mechanics Laboratory (DeeP MeLT), in collaboration with experimentalists from the Structural Engineering and Materials Laboratory (SEML): Dr. Lauren Stewart (Director), Mr. Jeremy Mitchell (Facilities Manager), Dr. Nan Gao (Research Engineer), Ms. Rebecca Tien (Graduate student). DEM simulations were conducted by Dr. Pei Wang, a former DeeP MeLT member and currently a research fellow at the University of Hong-Kong. The authors express their profound gratitude to all of these individuals for their assistance and support so far in this research project.
6

CHAPTER 1. INTRODUCTION AND BACKGROUND
1.1 Introduction During pre-stress transfer, the bottom portion of steel-reinforced concrete girders is
subjected to an important compression induced by the relaxation of tension in the bars. Girders tend to flex upwards and longitudinal cracks can develop along the horizontal axis. Reactions at the supports induce shear stress, which sometimes translates into additional diagonal cracks at the ends of the girders. During the subsequent lifespan of the girder, a variety of crack patterns can occur, including longitudinal (along the beam axis), transverse (perpendicular to the beam axis), and diagonal cracks. The National Cooperative Highway Research Program (NCHRP) recommendations are as follows (Tadros et al., 2010): Accept the girders if longitudinal cracks are less than 0.012 in. wide; Apply cementitious packing materials to cracks between 0.012 in. and 0.025 in. wide; Inject epoxy in cracks that are 0.025 in. to 0.050 in. wide; reject the girder (and replace it) if cracks are wider than 0.05 in. The Precast Concrete Institute (PCI) guide (2001) uses three defect categories: those that can be accepted without repair, those that can be repaired, and those that must be rejected. The PCI guide contains detailed trouble shooting and reparation procedures, depending on the location, orientation, length and width of the cracks. It specifies a minimum crack width for epoxy repair beginning at 0.006 inches.
Experimental studies showed that cracked concrete had lower mechanical strength than intact concrete, and that 80% of the lost concrete strength could be recovered by epoxy injection (Issa & Debs, 2007). It was also shown that epoxy injection can increase the fatigue resistance of cracked mortar, but that the mechanical performance of repaired concrete decreases as temperature increases (Shin et al., 2011). In general, repaired beams
7

exhibit a more brittle behavior than intact or corroded beams (Okada et al., 1988). According to the literature review completed in Phase I of this project, the model proposed by L. Bardella (Bardella, 2001) properly describes the nonlinear viscoelastic behavior of epoxy resins in the glassy state. But no model has ever been proposed to predict the behavior of cracked concrete repaired by epoxy. At present, there is a need to assess the mechanical integrity and sustainability of pre-stressed concrete beams during the entire life cycle of the built infrastructure, including crack propagation and reparation.
1.2 Prior Phases In Phase I of the project, we calibrated a Continuum Damage Mechanics model (the
Differential Stress-Induced Damage model, DSID) to simulate concrete crack propagation in the three directions of space. The DSID model is an energy-based model that predicts anisotropic crack propagation according to net tension and compression criteria (Xu and Arson, 2014; Xu et al., 2016). A static load test presented in (Tadros et al., 2010) was simulated with ABAQUS Finite Element (FE) software for four types of girders (Tennesse, Florida, Virginia, and Washington State designs), with and without pre-loading damage in the end zone. In order to study vertical cracks, horizontal cracks, diagonal cracks, and cracks on the bottom flange, damage was assigned to end zone elements in designated space directions prior to loading. The DSID model was assigned to the concrete elements and elastic bar elements were used to represent the pre-stressed steel reinforcements. The effect of initial cracks on deflection was minimal; increased displacements were only noted in the area close to the point of application of the clamping force. The most noticeable effects on deflection were obtained for vertical and diagonal cracks. It has to be noted that these results were obtained for smeared damage under a moderate static load. Then a
8

portion of typical bridge deck supported by pre-stressed girders was modeled with the Finite Element Method (FEM). Beam deflection and concrete damage induced by static loads were calculated with and without diffuse micro-cracking consecutive to pre-stress transfer. Initial damage configurations corresponded to typical crack patterns listed in the PCI bridge repair guidelines. In general, lower stress and higher deflections were found for higher initial damage.
The objective of Phase II is to simulate localized fracture propagation and epoxy reparation in concrete. Because experimental data on the mechanical properties of cracked concrete repaired by epoxy is scarce, the research objective of Phase III was to conduct experiments to measure the mechanical properties of cracked concrete repaired by epoxy in the laboratory, so as to calibrate and validate the numerical models proposed in Phase II. The research team developed and executed an experimental method to characterize epoxy-repaired concrete cylinders via uniaxial compression and splitting tension. The experiments included cutting plain concrete cylinders and reattaching them with epoxy. Uniaxial and splitting tension experiments were conducted on plain mortar cylindrical specimens in order to provide mechanical properties of the mortar as a separate material. Studies on specimen size effect were conducted in order to determine the most computationally efficient sample size for the Phase II models. Unreinforced concrete beam experiments were conducted on notched specimens. These tests included both plain concrete beams and beams that were repaired with an epoxy layer at a cut made directly above the notch. In all experiments involving the epoxy-repaired specimens, the crack initiated at the notch and then propagated through the concrete, away from the interface. The epoxy-repaired specimens exhibited slightly higher stiffness than the plain, notched
9

specimens. These experiments provided calibration data for the Phase II models, which allowed using the splitting test results for validation. Lastly, five reinforced concrete beam experiments were conducted on "as-built" and epoxy-repaired specimens in order to provide simple validation data at the system level. The epoxy-repaired specimens were loaded to induce cracking of various levels, unloaded, repaired with epoxy, and reloaded to failure. In the cases where the cracks were very small (< 0.006 in), the epoxy did not have a significant effect on the ultimate capacity, but the beam behaved in a more brittle fashion. Interestingly, in the experiments where the cracks were larger (> 0.006 in), the repaired beams exhibited a much higher capacity than the "as-built" ones. The failure mechanism and the ductility of the beam were both affected.
1.3 Phase II Objectives and Summary The objective of Phase II is to predict concrete mechanical recovery induced by epoxy injection, in order to check the girder acceptance criteria recommended by the NCHRP and by the PCI. Task 1: Numerical prediction of the mechanical performance of damaged concrete before and after epoxy injection. We combined a cohesive segment model with a continuum model of damage to predict the formation of macroscopic fractures in as-built and repaired concrete. We used the results of the mechanical tests performed in Phase III with plain and repaired concrete samples for model calibration and validation. This task also encompasses sensitivity analyses, aimed to understand the influence of the fraction of crack volume filled with epoxy on the mechanical behavior of repaired concrete.
10

Task 2: Assessment of the mechanical properties of repaired pre-stressed concrete girders. We implemented the numerical models proposed in Task 1 in numerical computation programs and to simulate the evolution of stress, deformation and damage in girders during pre-stressing and service loading. Four-point bending tests were simulated with and without pre-stress damage. Repair was modeled by replacing damaged elements by repaired elements in zones where the crack width exceeds a certain threshold after the loading phase. Four-point bending simulations were simulated for the repaired specimens. Lastly, the performance of as-built beams was compared to that of damaged beams and repaired beams, for four initial crack distributions described in the PCI guidelines, which allows assessment of the PCI recommendations.
1.4 Report Organization Chapter 2 of this report presents the models established and calibrated to understand
the mechanisms that drive fracture propagation in concrete. The Discrete Element Method (DEM) is used to understand the role of coarse aggregates in concrete strength, while the Concrete Damage Plasticity (CDP) is employed to predict fracture patterns at the laboratory scale, i.e. in plain concrete cylinders subjected to splitting and uniaxial compression tests and in plain concrete beams subjected to three-point bending tests.
Chapter 3 describes the multi-scale approach proposed to understand the behavior of plain concrete and reinforced concrete repaired by epoxy injection. To understand the nature of the bonding force between concrete and epoxy, we present a molecular-scale model of concrete/epoxy interface. Simulations of pullout tests under different environmental conditions and with different concrete mineral substrates shed light on the optimal conditions to use epoxy for concrete reparation. At the metric-scale, we propose a
11

method to model the reparation of concrete by epoxy injection with the FEM. The FEM model is enriched with cohesive zone (CZ) elements to represent the adhesive strength of cracks filled with epoxy. The constitutive parameters of the CDP and of the CZ models are calibrated to represent the properties of repaired concrete tested in Phase III. Parametric studies show that epoxy injection can be used to recover the stiffness and strength of structural members, at the condition that cracks that are wider than a certain width threshold are filled with epoxy. Experimental data gathered in Phase III is used for model calibration and validation.
Chapter 4 presents simulations of pre-stress, loading, repair and reloading in girders. The effect of pre-stress damage is analyzed. The performance of PSC girders with initial cracks is compared to PSC girders with repaired initial cracks. Results are discussed in light of the PCI recommendations for the types of initial crack distributions simulated.
Chapter 5 summarizes the work and provides some recommendations on the use of epoxy for concrete reparation.
Chapter 6 lists the references cited in this report.
12

CHAPTER 2. NUMERICAL SIMULATION OF FRACTURE PROPAGATION IN PLAIN CONCRETE
In this chapter, we first simulate crack propagation in concrete at the aggregate scale with the Discrete Element Method (DEM), to identify the types of defects that play the most critical role in fracture propagation (Section 2.1). Next, we simulate damage propagation and fracture pattern formation in concrete at the metric scale with the FEM, and we calibrate the CDP parameters against experimental results obtained in Phase III (Sections 2.2 and 2.3).
2.1 DEM Model: Role of coarse aggregates on concrete strength With a volume fraction of 40% to 50%, coarse aggregates play an important role in
concrete mechanical properties (Beushausen & Dittmer, 2015; Chi, Huang, Yang, & Chang, 2003). Zhou and collaborators found that the compressive strength of a cubic concrete sample is drastically reduced by weaker aggregates (Zhou, Lydon, & Barr, 1995). Similar results on compressive strength are reported in (Ahmad & Alghamdi, 2012; Aitcin & Mehta, 1990). Concrete strength is also greatly affected by the Interfacial Transition Zone (ITZ) between coarse aggregates and mortar. The ITZ is a region of 15 to 20 m in thickness, with a low content of cement particles and a high porosity, due to the wall effect (Scrivener, Crumbie, & Laugesen, 2004). The high local porosity and water content of the ITZ favors the deposition of calcium hydroxide, which, compared with calcium silicate hydrate (C-S-H), is endowed with weaker van der Waals forces and a more oriented microstructure, which provide less adhesion capacity. As a result, the ITZ is weaker in strength. The objective of this first study is to understand the micro-mechanisms that drive fracture
13

propagation in concrete and to assess the influence of the strength of the coarse aggregates as well as that of the ITZ properties on concrete strength.
2.1.1 Introduction to the DEM The DEM was proposed by Cundall and Strack in 1970s to perform research on the
behavior of granular assemblies (Cundall, 1971; Cundall & Strack, 1979). In this method, particulate interactions are determined by contact laws and particles' movements are governed by Newton's second law. We conducted the DEM simulations presented in this project with PF3D software. Rigid spheres and rigid walls are allowed to overlap at the contact. Contact surface areas are small compared to the size of the spherical elements. Spheres can be bonded together to form clusters of different sizes and shapes. As a result, the rheology of contacts is governed by three constitutive laws: a stiffness model; a slip model; and a bonding model. The stiffness model provides a relationship between a contact force and a relative displacement. The slip defines the maximum shear force within a contact so that slip movement can occur. The bonding model characterizes the rheology of the bonds, which can undertake forces and/or bending moments. Crack propagation is often represented by bond breakage, in which case, bonds disappear during the simulation.
In this project, we used the Hertz-Mindlin stiffness model in which the normal and shear stiffnesses increase with the sphere overlap (Itasca, 2008). The two governing parameters are the shear modulus G and Poisson's ratio . In the slip model, the maximum shear force at a contact is expressed as Fsmax=|Fni|, where is the friction coefficient and Fni is the normal contact force. We used the parallel bond model, in which bonds between spherical elements are represented as short beams of circular cross-section. Parallel bonds
14

can transmit both forces and moments between particles, and a detailed description of this model can be found in (Itasca, 2008).

2.1.2 Displacement softening contact bond model

The mechanical properties of the spherical elements that make a DEM model are the

normal stiffness kn, the normal/shear stiffness ratio kn/ks, the friction coefficient and the

mass density . Bonds between the spherical elements are brittle, which allows capturing

crack propagation in the DEM sample. Experimental results reported in numerous studies

indicate that, after the maximum load has been reached, concrete exhibits a gradual

decrease of loading capacity with increasing strain/displacement (Z. P. Bazant & Pfeiffer,

1987; Petersson, 1980). To capture this softening behavior, we propose a displacement-

softening contact model, inspired by the bond model proposed by Ma and Huang (Ma &

Huang, 2018). Figure 1 illustrates the normal and shear force-displacement curves

employed in the contact model, which is governed by the following five microparameters:

normal loading stiffness knl, maximum normal force Fnmax, normal softening stiffness kns,

bond stiffness ratio knl/ksl and maximum force ratio Fsmax/Fnmax. The maximum normal

force Fnmax is given by:

F!$ "#

=

(R%

+ R&)&' 4

=

A('

where R1 and R2 are the radii of the two particles in contact; st is the bond tensile strength;

A0 is the bonding area. The calculation is based on the assumption that the bonding area

between two spherical particles has a disc shape with a diameter equal to the average radius

of the two particles in contact. The maximum tensile force that a bond can carry is the

product of the bonding area by the tensile strength. The maximum force ratio Fsmax/Fnmax is
15

equal to the bond strength ratio t/ st. Therefore, two sets of micro-properties are used in

our DEM simulations, i.e. the micro-properties of the particles:

-k$

,

k$ k)

,

,

2

and the micro-properties of the displacement softening contact bond model:

k$*, ', k$), k$*/k)*, t/'. The relationships between the above micro-properties and the macroscopic properties of

the mortar and concrete DEM samples are discussed below.

Figure 1. Force-displacement relationships used in the displacement softening contact bond model
In our analysis, we followed the method proposed in (Potyondy & Cundall, 2004) and (Itasca, 2008) to define particle microparameters, as follows:
k$ = 4RE, k$/k) = E/G, = (1 + n)!"', where R is the particle radius; E and G are the Young's modulus and shear modulus of the material represented by the DEM sample, respectively; +,- is the mass density of the material and n is the porosity of the DEM specimen (around 35% in 3D simulations).
16

The normal loading stiffness knl is set equal to the average of the normal stiffnesses of the two bonded balls. The bond tensile strength is set equal to the tensile strength of the

material represented by the DEM sample. In our displacement softening contact model, the

bond shear strength is assumed to be higher than the tensile strength, so that shear stress

does not affect bond breakage. This is a realistic hypothesis for DEM models of bonded

and aggregated materials (Fakhimi & Villegas, 2007). The energy used to break a bond,

which is the area of triangle OAB in Figure 1, is equal to the product of the fracture energy,

Gf, by the bonding area, A0. The normal softening stiffness kns is expressed as:

k$)

=

(/!$ "#)%1$& &2'3(1$&4(/!$ "#)%

.

The remaining two contact bond parameters, knl/ksl and /- , control the shear

behavior of the material. In former studies (Ding & Zhang, 2014; Fakhimi & Villegas,

2007), a realistic ratio between the uniaxial compressive strength and the Brazilian tensile

strength (between 10 to 20 from most geomaterials) could be obtained with a high stiffness

ratio and strength ratio. Here, we calibrate knl/ksl and t/ st against experimental results. To

summarize, the required micro-properties can be expressed in the form:

E, G, , !"', ', G5, k$*/k)*, /'.

2.1.3 Mortar DEM model We generated cylinder specimens with the method proposed in (Potyondy & Cundall,
2004), which will be briefly introduced here. In a first step, an assembly of particles is generated in a material vessel. The number of particles generated depends on the volume of the vessel, an estimated porosity (35% in 3D) and the given average particle size. Initially, all particles are randomly located and may overlap. Then the radii of all particles

17

are modified uniformly so that an isotropic target stress, which is set to 1% of the uniaxial compressive strength, can be achieved at all boundaries. Large overlaps between particles are eliminated and force chains are generated. Particles with a coordination number (i.e., a number of contacts with neighboring particles) less than 3 are defined as ''floating'' particles in the sample. A densely packed sample can be generated by increasing the radius of these particles. Once the DEM sample is generated with the required porosity, the displacement softening contact bond model is assigned to all particle-to-particle contacts. Lastly, we remove the bounding walls and stress in the material is relaxed.
In our simulation, a cuboid mortar sample is first generated and then cut into the desired size and shape - in this case, 150100 cylinders for Brazilian Tests (BTs) and 100200 cylinders for Uniaxial Compression Tests (UCTs), like in the experiments. The minimum particle size was 1 mm and the size ratio between maximum particle size and the minimum particle size was 1.66. In each simulation, we had 120,000 to 130,000 particles and it usually took 12 hours on a workstation (i7-3770 3.4GHz CPU and 64 GB RAM) to simulate one test.
We obtain the mortar's Young's modulus from the empirical equation proposed in ACI 318 (Committee, 1999):
E = 4700?6 where f'c is the uniaxial compressive strength, reported in Table 1. We assume that the mortar's Poisson's ratio is 0.2. The shear modulus G is then given by E/(2(1+)). For the DEM friction coefficient, we use an empirical value of 0.5, as recommended in (Xu, Xu, Zhou, Du, & Hu, 2010). The tensile strength of the bonds, -, is set equal to that of the mortar specimen. Due to the size effects that occur in BTs, the tensile strength reported in
18

experiments is higher than the tensile strength of the material tested (Zdenek P Bazant, Kazemi, Hasegawa, & Mazars, 1991). We use the relationship proposed in (Rocco, Guinea, Planas, & Elices, 1999) to calculate mortar tensile strength from the BT strength and from the specimen diameter. For mortar, the fracture energy Gf ranges from 20 N/mm to 200 N/mm (Giaccio, Rocco, & Zerbino, 1993; Landis & Nagy, 2000). In our simulation, we use Gf = 80 N/mm, because this value provides the best fit with experimental results. With the DEM, it is challenging to simulate BTs and UCTs with a realistic ratio uniaxial compressive strength over Brazilian strength (UCS/BS). Different methods were used in previous studies, for example by generating angular particles to increase the interlocking forces, by increasing the initial compressive strength, by modifying the strength ratio, or by creating new contact models (Ding & Zhang, 2014; Fakhimi & Villegas, 2007). In our displacement softening contact model, the ratio of UCS/BS is mainly controlled by two parameters: knl/ksl and t/ st. In a Brazilian test, the failure mechanism is a combination of shear failure at the loading point and tensile failure at the center (D. Y. Li & Wong, 2013). A low knl/ksl decreases the shear stress at failure, thus yielding a lower BS. A large t/ st increases the shear strength measured in the uniaxial compressive test. Mortar parameters obtained after calibration are summarized in Table 1.
We conducted UCT and BT simulations twice to check that DEM randomization effects would not affect the calibration results. Simulation results are presented in Figures 2 to 5. The BTS and UCS are in agreement with the strengths measured experimentally, with an error of less than 5%. The ratio of UCS/BS, which is around 13.0 in our experiments, is properly reproduced in the DEM simulations. Figure 3 shows the location of broken bonds at several stages of the BTs. At 50% of the peak load, micro-cracks appear
19

near the loading platens, which is typical of a local shear mechanism close to the loading supports. At 90% of the peak load, the number of cracks near the loading platens increase. In addition, micro-cracks appear at the center, which reveals a tensile failure mechanism, induced by horizontal tensile stress. At the peak load, the micro-cracks that originated at the center propagate towards the loading platens. After the peak, the specimen rapidly breaks into two fragments.

Table 1. Mortar DEM model parameters

Parameter

Mortar

Young's modulus E (GPa):

36.8

Shear modulus G (GPa):

30.7

Friction coefficient (-):

0.5

Density mat (kg/m3):

2400

Tensile strength t (MPa):

3.85

Fracture energy Gf (N/m):

80

Stiffness ratio knl/ksl (-):

6.0

Strength ratio /t (-):

20

Aggregate 36.8 30.7 0.5 2400 3.5 72.7 6.0 20

Figure 2. Stress-strain curves of mortar in the Brazilian tests simulated with the DEM
20

Figure 3. Micro-cracks in mortar at different stages of the Brazilian tests
In the UCT, we observe randomly distributed micro-cracks at 50% of the peak load. No obvious shear band or fracture can be identified at this point. When the load reaches 90% of the peak load, we can observe more micro-cracks, preferentially oriented along the directions of two shear failure planes (shear bands). At this loading stage, the stiffness of the material starts to decrease due to material softening induced by micro-crack propagation. As we continue loading from 90% to 100% of the peak load, the number of broken bonds in the shear bands increases rapidly. Figure 5 clearly shows a fracture pattern characterized by more than one oblique plane. This failure mode is one of the four most common failure modes reported by ASTM codes (ASTM), along with diagonal fractures,
21

columnar vertical cracks and the combination of conic and vertical cracks. The observed failure mechanism is typical of fixed ends boundary conditions, in which strong confinement is used at loading platens to prevent any lateral displacement of the specimen. The fixed ends boundary condition is used in both our experiments and simulations.
Figure 4. Stress-strain curves of mortar in the uniaxial compression tests simulated with the DEM
22

Figure 5. Micro-cracks in mortar at different stages of the uniaxial compression tests
23

2.1.4 Concrete DEM model In order to generate realistic aggregate shapes, we scanned 20 representative
aggregates with an X-ray scanner. Aggregate shapes were characterized by 20,000 to 30,000 points, which provided detailed information about particle shapes and surface textures. Figure 6 shows the original scan data and the corresponding DEM shapes obtained for four aggregates.
Figure 6. Generation of aggregate shapes in the DEM model
The DEM concrete model is generated as follows: 1. Create the mortar sample using the method presented in Section 2.1.3. Calculate the number, sizes and positions of the aggregates. The number of coarse aggregates in each fraction size, ni is given by n7 = -(&8 - %8)/8 where Vt is the total volume of the specimen; is the volume fraction of the coarse aggregates; P1i and P2i are the percentage of aggregates finer than d1i
24

and d2i in the PSD curve, respectively; Vi is the average volume of an aggregate in this size fraction. Within each size fraction, the size of an aggregate is assumed to follow a uniform distribution. A MATLAB code was written to generate non-overlapping spheres within the sample boundary, with the given number of aggregates and the given aggregate sizes. Each sphere generated by the MATLAB code corresponds to an aggregate of specific volume. 2. Replace the spheres by realistic aggregate shapes. We replace each sphere with a realistic shape randomly chosen from the set of 20 shapes obtained by X-ray scanning. After the substitution, some aggregates may overlap, which slightly decreases the total volume fraction of the aggregates. In addition, some aggregate fractions may be outside of the boundary of the specimen. 3. Identify the spherical elements and the contacts that belong to the aggregates. Another MATLAB algorithm is created to loop through all spherical elements and contacts, and check whether they belong to an aggregate or not. Then the aggregate volume fraction is calculated as the ratio of the number of spheres in the aggregates by the total number of spheres in the specimen. We compare this volume fraction with the aggregate volume fraction found experimentally, update the parameter accordingly, and repeat steps 2 to 4 if the difference is larger than 1%. 4. Change the properties of the spherical elements and bonds within the aggregate contours and within the ITZ. To account for the local increase of
25

porosity in the ITZ, we removed a certain portion of displacement softening contact bonds at the aggregate/mortar interface. 5. Generate voids. Another important difference between mortar and concrete is the void ratio. An estimated void ratio of 4% is used in our concrete model, which is a reasonable estimation according to images of X-ray Computed Tomographic images obtained in (Gallucci, Scrivener, Groso, Stampanoni, & Margaritondo, 2007). Figures 7(a) and (b) show the DEM concrete samples used to simulate the concrete BTs and UCTs, respectively.
Figure 7. DEM concrete specimens used to simulate (a) Brazilian tests and (b) Uniaxial compression tests
26

The coarse aggregates used in the concrete tested in our experiments are made of granite from Norcross, GA. In some concretes, coarse aggregates are found to be much stronger than mortar and have a much lower probability to break (Skarzynski & Tejchman, 2016). Contrary to those observations, in both the BTs and UCTs we conducted, a large number of coarse aggregates were crushed and fresh aggregate surfaces were exposed, which indicates that aggregates had a lower tensile strength and a lower compressive strength than mortar. Not surprisingly, our calibration simulations yielded a lower bond strength - for coarse aggregates than for mortar. In addition, we assume that the ultimate bond displacement at failure is the same as that of mortar, so that the fracture energy of the aggregates, 9, is proportional to the aggregate bond strength. The ITZ also contributes to the relative lower strength of concrete compared to that of mortar (Scrivener et al., 2004). The existence of coarse aggregates leads to a local increase in porosity and thus to a smaller bond area between aggregates and mortar. To account for the effect of the ITZ, we removed 20% of the displacement softening contact bonds between aggregate particles and mortar particles. This value of 20% was obtained by trial and error (the effect of this ratio on concrete strength is discussed in the following).
We conducted the BT and UCT simulations three times, to check if DEM randomization would have any effects on the results. Simulation results are shown in Figures 8 to 11. The coarse aggregate properties obtained after calibration are summarized in Table 1. Concrete stress-strain curves are similar to those of mortar. However, because the coarse aggregates increase the heterogeneity of the specimens, a variability of 16.0% (respectively 7.0%) is observed for the BTS (respectively UCS) among the three simulations (Figures 8 and 10). The sequence of bond breakage events in concrete is similar
27

to that in mortar. In the BT, micro-cracks first appear at the loading platens but the failure is caused by tensile micro-cracks that originate at the center of the specimen and then propagate towards the platens. The failure plane after the peak load is not a straight line, because of the specimen heterogeneity induced by the presence of aggregates. Like in mortar, multiple shear failure planes are observed, see Figure 11. Due to the lower loading capacity of concrete, the number of micro-cracks is less than in mortar. In other words, fractures form at an earlier stage, passing through the weak aggregates and the ITZ. A detailed analysis of the influence of the relative fraction of inactive bonds at the ITZ is presented in the next section. To sum up, the existence of coarse aggregates greatly affects both the BTS and the UCS, but the presence of aggregates does not influence the failure mechanisms involved.
Figure 8. Stress-strain curves of concrete in the Brazilian tests simulated with the DEM
28

Figure 9. Micro-cracks in concrete at different stages of the Brazilian tests
Figure 10. Stress-strain curves of concrete in the uniaxial compression tests simulated with the DEM
29

Figure 11. Micro-cracks in concrete at different stages of the uniaxial compression tests
30

2.1.5 Sensitivity analyses
In Section 2.1.4, we calibrated our DEM model to reproduce the behavior of concrete, which exhibited lower strength than mortar. Comparing a mortar and a concrete of same water/cement ratio and same aggregate volume fraction (fine aggregates in mortar and coarse aggregates in concrete), it is found that mortar has a higher strength than concrete. This observation led some authors to think that the ITZ acts as a plane of weakness in concrete (Scrivener et al., 2004). The ITZ is a layer of around 30 m that can only be observed through scanning electron microscopy (SEM); it is hence challenging to monitor the ITZ during concrete loading. In addition, it was noted that concrete strength decreases as the strength of coarse aggregates decreases (Chi et al., 2003), but most studies focused on strong aggregates materials. Despite qualitative observations, quantitative analyses of the effect of the ITZ and aggregate strength on concrete strength are lacking. Therefore, in this section we study the sensitivity of concrete strength to the ITZ and to aggregate strength by simulating the BT and the UCT with our DEM model. We generate the DEM specimens following the method explained in Section 2.1.4 and we use the calibrated parameters listed in Table 1.
2.1.5.1 Effect of aggregate shape Previous research shows that concrete containing angular crushed rock is stronger
than concrete containing smooth gravel. In our model, aggregates are modeled with weaker bonds and softer particle elements, and shapes are determined by X-ray scanning. In order to evaluate the effect of aggregate shape on concrete strength, we simulate the response of a DEM mortar sample with the same fractions of broken bonds (ITZ bonds), weaker bonds (aggregate inter-particle bonds) and softer particles (aggregate particles) as those of the
31

calibrated DEM concrete samples, but with a random distribution of broken bonds, weak bonds and soft particles. The effect of aggregate shape on concrete strength is shown in Figs. 12 and 13. When aggregate shape is ignored, the specimen can be regarded as a homogeneous material with uniformly distributed micro cracks and flaws. By contrast, in the concrete model that accounts for scanned aggregate shapes, the cracks and flaws are localized and higher stress concentrations occur. As a result, strength ranks as follows: mortar > concrete with uniformly distributed flaws (aggregate shape ignored) > concrete with scanned aggregate shapes. Simulation results show that the strength of concrete with scanned aggregate shapes is around 20% lower than of concrete with aggregate shapes ignored. This result indicates that concrete is stronger when aggregate/mortar interfaces are uniformly distributed, which is more likely to occur when the aggregate size is smaller.
5

4

Stress (MPa)

3

2

1

Mortar

Concrete with scanned aggregate shapes

Concrete with aggregate shapes ignored

0

0

1

2

3

4

5

Strain

10-3

Figure 12. Effect of aggregate shape on concrete strength in Brazilian tests.

32

70

60

Stress (MPa)

50

40

30

20

Mortar

10

Concrete with scanned aggregate shapes

Concrete with aggregate shapes ignored

0

0

1

2

3

4

5

Strain

10-3

Figure 13. Effect of aggregate shape on concrete strength in uniaxial compression tests.

2.1.5.2 Effect of the Interfacial Transition Zone (ITZ) on concrete strength In our model, the weakening effect of the ITZ is accounted for by deleting a fraction
of the aggregate/mortar bonds. We define the contact ratio, , as the number of displacement-softening contact bonds used at aggregate/mortar interfaces over the total number of interfacial contacts. The model was calibrated for =80%. We now vary between 50% and 100% to understand the influence of mortar/aggregate adhesive surface on concrete BTS and UCS. Results are presented in Figures 14 and 15, respectively. Note that for each contact ratio we conducted the BT and UCT simulations at least three times. Only the average results are shown here. The ITZ greatly affects concrete strength: the BTS and the UCS respectively increase by 227.5% and 222.0% when varies from 50% to 100%. The stress/strain curves are similar in all simulations: At first, stress increases with strain but the stiffness gradually decreases due to bond softening and breakage; Then

33

stress reaches a peak value and the specimen starts to fail; As the specimen is further compressed, stress decreases rapidly and numerous bonds break. Figure 16 and 17 show that the strength-contact ratio relationship is linear for both the BT and the UCT. The higher , the higher the strength. Contacts without displacement softening contact bonds can be viewed as internal micro-cracks or micro-flaws, similar to the ITZ in actual concrete. Micro-crack interaction and coalescence occur in specimens with low contact ratio, which drastically lowers concrete strength.
Figure 14. Stress-strain curves of concrete with different interface contact ratios in the Brazilian tests
34

Figure 15. Stress-strain curves of concrete with different interface contact ratios in the uniaxial compression tests
Figure 16. A quasi-linear relationship between BTS and contact ratio
35

Figure 17. A quasi-linear relationship between UCS and contact ratio
2.1.5.3 Effect of aggregate tensile strength on concrete strength In the concrete tested in this study, the volume fraction of coarse aggregates amounts
to nearly 45%. Aggregates are thus expected to play an important role in concrete mechanical behavior. The calibrated tensile strength of the coarse aggregates is 3.5 MPa, which is lower than the mortar BTS. During the BTs and UCTs, we note indeed that concrete strength is lower than that of mortar, and we observe a large number of crushed aggregates. Materials commonly used for coarse aggregates include quartzite, limestone, marble, low strength granite and high strength granite. Corresponding aggregate tensile strengths span from 3.0 MPa up to 15.0 MPa. In order to understand the effect of aggregate tensile strength on concrete BTS and UCS, concrete specimens with different aggregate strengths are generated and subjected to BTs and UCTs. As mentioned in Section 2.1.4, we assume that the ultimate bond displacement at failure is the same for all aggregates, so that the change of aggregate strength dictates the change of fracture energy, 9. BT and UCT
36

simulation results are shown in Figures 18 and 19, respectively. In general, concrete strength increases with aggregate strength, except in the simulation of the UCT with coarse aggregates with a tensile strength of 3.5 MPa. The exact reason for this unknown; this discrepancy may stem from the variability of concrete specimens generated with the DEM. Figures 20 and 21 show that the relationship between aggregate tensile strength and concrete strength is quasi-linear both on the BT and in the UCT. It is also noticeable that when the aggregate strength is increased by more than 300%, concrete BTS and UCS increase by less 68%. In other words, both the ITZ and aggregate tensile strength influence concrete strength, but concrete strength is most sensitive to the contact ratio in the ITZ.
Figure 18. Stress-strain curves of concrete with different aggregate tensile strengths in the Brazilian tests
37

Figure 19. Stress-strain curves of concrete with different aggregate tensile strengths during the uniaxial compression tests
Figure 20. A quasi-linear relationship between BTS and coarse aggregate tensile strength
38

Figure 21. A quasi-linear relationship between UCS and coarse aggregate tensile strength
2.1.6 Summary and conclusions We constructed DEM models of mortar and concrete by generating samples of
bonded particles. First, the mechanical parameters of the spherical elements and of the displacement softening contact bonds were calibrated against mortar splitting and uniaxial compression tests. Then, coarse aggregates were generated from point clouds obtained by scanning actual aggregates used in concrete. The aggregate bond parameters were calibrated against concrete splitting and uniaxial compression tests, with the calibrated mortar parameters. The proposed DEM models capture the mechanical behaviors of mortar and concrete, with a realistic ratio of uniaxial compression strength / splitting strength. We studied the influence of the tensile strength of coarse aggregates, the spatial distribution of aggregate/mortar interfaces and the adhesive area of the ITZ on concrete strength. The main conclusions are the following:
39

1. Concrete splitting and uniaxial compression strengths increase with aggregate tensile strength and with the ITZ contact ratio .
2. At same aggregate volume fraction, concrete strength increases when aggregate size decreases.
3. Mortar splitting and uniaxial compression strengths were about twice as much as concrete splitting and uniaxial compression strengths, respectively. Other authors had also noticed that concrete could have lower strength than mortar, and attributed this phenomenon to the weak surfaces in the ITZ. Our study confirms that the ITZ is the primary cause of low concrete strength. The aggregate volume fraction is usually around 30% or more. We found that multiplying the aggregate tensile strength by four could increase concrete strength by 2/3.
4. The proposed displacement softening contact bond model can capture mortar and concrete failure mechanisms in both the splitting test and the uniaxial compression test. In the splitting test, micro-cracks first appear at the vicinity of the loading platens. But failure is induced by fractures that initiate at the center of the specimen and that coalesce with the cracks located close to the platens. Specimens subjected to uniaxial compression tests exhibit multiple shear failure planes -- usually, two symmetric shear bands. A realistic ratio of uniaxial compression strength / splitting strength is obtained. Micro-parameters used in the proposed bond model are closely related to the macro-scale material properties, which facilitates calibration.
5. The ITZ is accounted for by deleting bonds at the mortar/aggregate interface. The ITZ can be seen as a distribution of interfacial micro-cracks, like in actual concrete. Both
40

in the splitting test and in the uniaxial compression test, a linear relationship exists between the contact ratio in the ITZ and concrete strength, for 0.5 < < 1.0.
6. Concrete strength is linearly related to aggregate tensile strength. Weak aggregates decrease concrete strength. However, the influence of aggregate tensile strength on concrete strength is much lower than that of the ITZ.
2.2 CDP Model: Fracture patterns in plain concrete 2.2.1 Introduction to the CDP model
We used the Concrete Damage Plasticity (CDP) model built in ABAQUS standard software to simulate crack propagation in plain, reinforced and pre-stressed steelreinforced concrete with the FEM. The CDP model was developed based on previous research conducted by Lubliner and collaborators (1989) and Lee and Fenves (1998) for modeling concrete and quasi-brittle materials. The CDP model couples the concepts of damage elasticity with tensile and compressive plasticity to account for the inelastic behavior of concrete, as shown in Figure 22. In the plots presented, !cpl and !tpl refer to the compressive and tensile plastic strains, respectively. The stress-strain relations are given by:
c = (1 - dc )E0 (c - !cpl ) t = (1 - dt )E0 (t - !tpl )
where dc and dt are the compression and tension damage variables, respectively. E0 is the elastic modulus of the material in the non-damaged state. The damage variables can be calculated as:
41

dc

=

c c-max

dt

=

t t-max

Figure 22. Principle of the CDP model built in ABAQUS (a) compressive plasticity (b) tension plasticity

Tensile and compressive inelastic deformation and damage parameters are

determined from the splitting tension and uniaxial compression tests performed in Phase

III, respectively. The tensile and compressive strengths of the 28-day concrete were found

to be 2.46 MPa and 35.54 MPa in Phase III experiments, respectively. The stress-strain

response in compression is assumed to obey the following constitutive relationships

(Mander et al, 1988):

6

=





: 66



- 1 + ;

= 6 66

42

66 = 6< J1 + 5 L::66<6 - 1MN

= 6 6 - =>6

=>6

=



: 66

66

68? = 6 - 6>@,=-86

Where

6

is

the

uniaxial

compressive

stress,

6

denotes

corresponding

strain,



: 66

is

the

compressive strength, and 6< is a set reference value of 0.02. The compressive stress-

strain curve can be obtained by calculating the mechanical state at different points, as

explained in Table 2.

Stress [MPa] 26.4194 27.8352 29.0449 30.054 30.8712 31.508 31.9776 32.2945 32.4736 32.53 32.4784 32.3328 32.1064 31.8113 31.4586 31.0582 30.6191 30.1491 29.655

Table 2. CDP input stress-strain curve in compression

Inelastic Damage

Stress

Inelastic

strain parameter [MPa]

strain

0

0

29.1429 0.00278895

0.000335891

0

23.8588 0.00428895

0.000422261

0

19.3844 0.00578895

0.000519185

0

15.9681 0.00728895

0.000626199

0

13.3899 0.00878895

0.000742704

0

11.4195

0.0102889

0.000868004

0

9.88521

0.0117889

0.00100134

0

8.66756

0.0132889

0.00114192

0

7.68392

0.0147889

0.00128895

0

6.8766

0.0162889

0.00143895 0.00158637

6.2046

0.0177889

0.00158895 0.00606182 5.63822

0.0192889

0.00173895 0.0130219

5.15558

0.0207889

0.00188895 0.0220933

4.74024

0.0222889

0.00203895 0.0329358

4.37969

0.0237889

0.00218895 0.045243

4.06425

0.0252889

0.00233895 0.0587418

3.78631

0.0267889

0.00248895 0.0731909

3.53986

0.0282889

0.00263895 0.0883789

Damage parameter
0.104122 0.266561 0.404106 0.509127 0.588383 0.648956 0.69612 0.733552 0.76379 0.788607 0.809265 0.826676 0.841513 0.854281 0.865364 0.875062 0.883606 0.891182

43

The tensile behavior of concrete is represented by a bilinear softening curve (see

Figure 22 above), which shows good agreement with the load-displacement curve obtained

experimentally in previous studies (Wittmann et al 1988, Guinea et al 1994, Park et al

2015). To find the tensile strength of the concrete material, we calculated the average

indirect tensile strength measured from eight cylinders samples with dimension of 150mm

diameter and 300mm height in Phase III. The average value of indirect tensile strength

(8-=) was 2.46MPa. To determine real-tensile strength (;-=) which controls the microcrack and damage initiation in plain concrete, we introduce the following relationship

(Rocco et al. 1999):

8-=

=

;-= -17.57 + 66.05

6 @

+

1.01305;-=

@

=

6 6 2;-=

where @ is a characteristic length, 6 is the diameter of the cylinder specimen and 6 is the

intercept of the initial tangent to the bilinear softening curve, as proposed by Elices et al.

(2002). The real-tensile strength is size independent (Coronado and Lopez 2007). Based

on the experimental results of Phase III, it is found be 2.26 MPa. The fracture energy is

calculated based on the empirical equation proposed by Bazant and Becq-Giraudon (2002)

and then modified by Martin et al. (2007) as follows:

9

=

0.0025(

W

6

(.B
Y

0.058

W1

+

+,C Y(.BD 1.94

4(.%E \^


where (is an aggregate shape factor (angular aggregate:1.12 and rounded aggregate:1);

6

is

the

compressive

strength;

+,C is

the

maximum

aggregate

size;

F 6

is

the

water-cement

ratio. We found that 9 is 0.105N/mm, which is within the range of values expected for

44

concrete fracture energy, which is typically between 0.084 and 0.13N/mm (Martin et al. 2007). We choose a stress ratio of 1.2 at the kink point () based on previous studies (Bazant and Becq-Giraudon, 2002). The real-tensile strength, fracture energy, and stress ratio at the kink point are used to calculate the damage parameters and to relate stresses and strains, as shown in Table 3.

Table 3. CDP input stress-strain curve in tension

Stress [MPa]

Damage Cracking displacement
parameter

2.28

0

0

1.883931429

0.01

0.173714286

1.487862857

0.02

0.347428571

1.091794286

0.03

0.521142857

0.695725714

0.04

0.694857143

0.456

0.046052632

0.8

0.371965714

0.08

0.836857143

0.272948571

0.12

0.880285714

0.173931429

0.16

0.923714286

0.099668571

0.19

0.956285714

0.074914286

0.2

0.967142857

We choose a standard Poisson's ratio of 0.2 for concrete. The Young's modulus is calculated based on the model of Mander et al. (1988):

6 = 5,000? 6( ( )

in which 6( is the compressive strength. Additional plasticity parameters include the

dilation angle, eccentricity, the ratio 9)*
9+*

and the coefficient Kc respectively chosen as 30

0.1, 1.12 and 0.75, based on previous studies (Jankowiak 2005; Labibzadeh et al, 2017).

Table 4 summarizes the concrete material parameters used in the CDP model. Simulations

presented in the next Section aim to verify the CDP model.

45

Table 4. Concrete material properties in the CDP model

Elastic

modulus

Poisson's

Dilation Eccentricity

G<

(GPa)

ratio

angle

6<

Kc

26.608

0.2

30

0.1

1.12

0.75

2.3 Verification of the CDP model against plain concrete experiments 2.3.1 Uniaxial compression tests on plain concrete cylinders
In order to verify the CDP in compression, uniaxial compression tests conducted in Phase III on plain concrete are simulated. Concrete cylinders were 12-inch height and had a 6-inch diameter. The experimental set-up and the corresponding numerical model are presented in Figure 23.

a

b

Figure 23. (a) Experimental setup of the uniaxial compression test (b) FEM model employed to simulate the uniaxial compression tests.

46

The rigid platens at the top and bottom of the concrete sample were modeled explicitly with steel Finite Elements. The steel/concrete interface was modeled with a hard contact (non-penetrability condition) with a friction coefficient of 0.1. Fixed displacements were imposed to the centroid of the bottom platen. A vertical displacement was imposed to the centroid of the top platen, until the ultimate axial compressive displacement was reached. A linear hexahedral element of type C3D8 with a mesh size of 5mm was used in the simulations. Figure 24 presents the stress-strain curve obtained by simulating the uniaxial compression test (black dots) and the stress-strain curve that can be obtained empirically from the experimental results, by using the constitutive relationships presented in the previous section (blue line). Note that the numerical stress-strain curve was obtained by extracting the force and displacement at the node located at the center of the top rigid platen. The red dashed line marks the compressive strength measured experimentally (32.54 MPa). Both the numerical and the empirical stress-strain curves match the experimental strength. The numerical and empirical stress-strain curves match.
Figure 24. Stress-strain curve obtained by simulating the uniaxial compression tests with the CDP model (black dots) in comparison with the input constitutive law (blue line) and
with the experimental strength (red dashed line).
47

Figure 25 shows that the distribution of compression damage forms clear shear bands and
that the concrete cylinder ultimately exhibits a columnar failure, like in the Phase III
experiments.
DAMAGEC (Avg: 75%)
+2.056e+00 +1.000e+00 +9.167e01 +8.333e01 +7.500e01 +6.667e01 +5.833e01 +5.000e01 +4.167e01 +3.333e01 +2.500e01 +1.667e01 +8.333e02 +0.000e+00

ODB: norepairerate01.odb Abaqus/Explicit 6.121 Fri Jan 31 17:32:48 Eastern Standard Time 2020

Step: Step1

Y

X

Increment 113695: Step Time = 0.3000

Primary Var: DAMAGEC

Deformed Var: U Deformation Scale Factor: +1.000e+00

Z

Figure 25. Damage distribution during the uniaxial compression test

2.3.2 Notched three-point bending test without reinforcement We simulate the notched three-point bending test performed during Phase III to
verify that the fracture energy used in the CDP model correctly represents the behavior of concrete in tension. The experimental set-up and the numerical model are explained in Figure 26. The beam was 1219 mm long, 250 mm high and 150 mm deep. The vertical load (displacement controlled) was applied to a rigid steel plate on the center-top of the
48

beam. The mesh size around the pre-existing notch is refined to avoid mesh sensitivity. We use C3D4 type elements in a 3D model.
Figure 26. Three-point bending tests on notched plain concrete beams. Top: Experimental set-up. Bottom: numerical model.
The distribution of damage in the beam is represented in Figure 27. As expected, tensile damage propagates vertically ahead of the notch. In Figure 28, Force-CMOD curves obtained in the laboratory are compared to Force-CMOD curves obtained for different kink point stress ratios and tensile strengths, showing that a satisfactory match can be obtained for the constitutive parameters chosen in this study (stress ratio of 1.2 and tensile strength of 2.2 MPa). The peak load in the experiment was 12.5 kN, and that predicted in the simulation is 14 kN, which confirms the suitability of the CDP model for simulating damage propagation in concrete.
49

Figure 27. Distribution of tensile damage during the three-point bending test simulations conducted with the CDP model under a vertical displacement of 0.288 mm (left) and 3.08
mm (right).
Figure 28. Numerical and experimental force-CMOD curves obtained during the threepoint bending tests. Left: Sensitivity of the model to the stress ratio at the kink point. Right: Sensitivity of the model to the tensile strength.
2.3.3 Three-point bending test on reinforced concrete beams The reinforced concrete (RC) beams tested in Phase III were fabricated with the
same dimensions as the non-reinforced beams, without any initial notch. Two #6 steel
50

rebars with a yield strength of 420 MPa were placed longitudinally as shown in Figure 29. We used a string potentiometer to measure the vertical displacement at the bottom center of the beams. Additionally, we speckled a longitudinal face of the beams and we measured the displacement field on that face by using the Digital Image Correlation (DIC) software developed by Blaber and collaborators [26]. Two RC beams were tested "as built", to investigate the force-displacement behavior without epoxy injection.
1200 mm
250 mm

1000 mm

150 mm

Figure 29. Geometry of the reinforced concrete beam. The rebars were placed at a distance of 45 mm from the bottom of the beam and the spacing was 60 mm.

We represent the steel rebars with linear elastic truss elements, with a Young's modulus of 200 GPa and a Poisson's ratio of 0.2. The distribution of damage is illustrated in Figure 30. In the early stages of the loading, several vertical damaged zones develop at the bottom center of the beam, due to tension. Close to failure, large diagonal cracks emerge from the supports, due to shear. These trends reflect the observations made in the experiments, in which multiple vertical cracks were followed by major shear cracks emerging from the support and going towards the center of the beam. In addition, the predicted force-displacement curve shows a good agreement with the experimental results obtained for both specimens tested in the laboratory, as can be seen from Figure 31.

51

Figure 30. Damage evolution during the three-point bending test conducted on reinforced concrete beams. Top: Numerical results (under displacements of 1.9 mm and
4mm). Bottom: Photograph of the experiment showing a specimen at failure.

150 RC concrete #1 RC concrete #2 FE simulation
100

Load [kN]

50

0

0

2

4

6

8

10

Displacement [mm]

Figure 31. Force displacement curves obtained numerically and experimentally for the three-point bending test conducted on reinforced concrete beams

52

2.3.4 Summary and conclusions We simulated laboratory experiments conducted in Phase III on GDOT concrete by
employing the CDP model in ABAQUS FEM software. The mechanical parameters of the CDP model were calibrated against splitting, uniaxial compression and three-point bending tests conducted on non-reinforced concrete. We then verified that the calibrated CDP model could be used to simulate the three-point bending tests conducted on reinforced concrete beams. We showed that the choice of the stress ratio at the kink point had little influence on the load-displacement curves obtained numerically. The CDP captures the average mechanical behavior of the specimens through the force-displacement curves, as well as the fracture patterns observed experimentally.
53

CHAPTER 3. MULTI-SCALE STUDY OF PLAIN CONCRETE REPAIRED BY EPOXY DEPOSITION
This chapter presents a study of concrete/epoxy interfaces at molecular and metric scales. We first demonstrate the feasibility of using Molecular Dynamics (MD) to simulate the tensile behavior of silica/epoxy interfaces, by comparing our simulation results to data published on the most common epoxy studied in the literature (DGEBA with DETA). We then simulate pullout tests on interfaces between the epoxy used by GDOT (High Molecular Weight Metacrylate, HMWM) and various concrete minerals, under several conditions of temperature and humidity. MD results shed light on the nature of the forces responsible for interface strength and indicate the optimal environmental conditions for repairing concrete by epoxy deposition. At the metric scale, we use the CDP model calibrated in Chapter 2 to simulate laboratory tests conducted on repaired concrete in Phase III. The FEM CDP model is enriched with Cohesive Zone (CZ) elements that represent the cracks filled with epoxy. Results allow quantification of stiffness, strength and load capacity of specimens before and after reparation.
3.1 MD model of concrete/epoxy interface 3.1.1 Introduction to MD
Interfaces between organic and inorganic materials are important research topics in chemistry, material sciences, mechanical engineering and civil engineering. Finite Element Methods and Finite Difference Methods were used on the basis of classical continuum mechanics theories. However, these approaches fail to capture structural interlocking and
54

intermolecular non-bonded forces at bi-material interfaces, which can only be modeled at the molecular scale. According to the early works of Sasse and Fiebrich (Sasse & Fiebrich, 1983), the strength of biomaterial interfaces is indeed caused by mechanical interlocking induced by the percolation of one material into the voids of the other (Figure 32) and by Van der Waals forces.
Figure 32. Examples of interface interlocking (Sasse & Fiebrich, 1983)
MD simulations, first introduced by Alder and collaborators (Alder & Wainwright, 1959), allow solving for the inter-molecular interaction forces that exist between organic and inorganic materials. MD models represent bi-material interfaces at the atomic scale, i.e. the length and time scales used in MD differ by several orders of magnitude from those of the laboratory experiments used to determine interface mechanical properties. Several up-scaling techniques were proposed in order to use MD simulation results for predicting macroscopic cohesion. For instance, Lau and collaborators calculated the intrinsic strength and fracture energy of an individual epoxy chain and calculated the macroscopic interface fracture energy by calibrating an epoxy chain spacing parameter (Lau, Bykztrk, & Buehler, 2012). This interface fracture energy is then used in a cohesive segment model. This rubber-like elasticity model is illustrated in Figure 33. The main drawback of this
55

model is its dependence to the discrete grid size, which controls the spacing (s). In this project, we will present two separate studies at molecular and metric scales, and draw complementary conclusions from each scale of observation.
Figure 33. Multiscale analysis of epoxy-silica interface: "rubber-like model" (Lau et al., 2012)
The steps to build an MD model are explained in the flow chart shown in Figure 34. In the initialization step, the molecular structure of the materials is optimized, the boundary conditions are applied and the potential energy function is defined. In the simulation step, Newton's second law is used to calculate how forces will affect the motions of atoms, as follows:
() () =

56

= () = - () Where F is the force exerted on an atom at time t, m is the mass of the atom, a is the atom's acceleration, and x is the atom's position. U is the potential energy function. Verlet algorithm, based on the Finite Difference Method, is often used in MD to integrate the equations of motion.
Figure 34. MD simulation flow chart (Lau et al., 2012)
Verlet algorithm is given by: 8(( + ) = 28(() - 8(( - ) + 8(()()&
where xi(t0) is the position vector at t0, xi(t0-Dt) is the position vector at t0-Dt and a(t0) is the acceleration vector at time t0. Atomic positions and velocities are updated iteratively, based on the knowledge of the initial positions and velocities, and of the potential energy function. The latter depends on the forces experienced by any atom given the positions of the other atoms. The choice of the force field function plays an important role in the
57

calculation of the potential energy of the molecular system. The force field is typically

derived from experimental datasets and approximations. Classically, the potential energy

function is composed of a bonded force field potential and of non-bonded force field

potential, as follows:

-<-,@ = G<?H>H + ?<?G<?H>H

G<?H>H = l G( - ()& + l I( - ()& + l ()&

G

I



?<?G<?H>H

=

l 8K
8,K

tL8K

%&
M

8K

-

L8K 8K

M
M

u

,

8K



6

The superposition of bond stretch, bond angle, and dihedral angle torsion terms in the

second equation represents the covalent bond energy, in which G , I , and are constants, and b0, q0 and are the equilibrium bond length, bond angle, and dihedral angle. rij is the distance between the ith and jth atoms. The non-bonded energy term is the Van der Waals energy (Lennard-Jones potential). The non-bonded force energy potential is the

most computationally expensive in MD models.

Before conducting any MD simulation, not only the force field function has to be

defined and validated against experimental data, but also, the energy of the system has to

be minimized, and the boundary conditions have to be properly defined. The optimization

of the geometry of the system, i.e. the system's energy minimization, is achieved by solving

a partial differential equation, like in continuum mechanics. Atomic coordinates are

calculated to ensure that the potential energy reaches a local minimum. The derivative of

the energy relative to the atoms' positions \OP^ is set to zero, and the second derivative of
OC

the energy potential function relative to the atoms' positions W OOP Y is constrained to be
OC,OC-

58

always positive. Additionally, it is required to ensure consistency between the ensemble that is used to calculate the admissible microscopic configuration and that that represents the macroscopic state of the materials. The three ensembles commonly used in MD are the micro-canonical ensemble (NVE), the canonical ensemble (NVT) and Gibb's ensemble (NPT). Table 5 shows the possible ensembles: each ensemble comprises one and only one field variable from each row.

Table 5. Possible Ensembles in MD Particle number (N)

Volume (V)

Pressure (P)

Energy (E)

Temperature (T)

MD was widely used for modeling epoxy materials. The simulation of the cross-linking of epoxy with an epoxy hardener is one of the most challenging tasks in MD simulations. Bandyopadhyay and collaborators constructed a MD model for a particular epoxy polymer with a cross linking agent in order to understand the influence of the stoichiometric mixture of the epoxy (EPON-862) and of the hardener (DETDA) as well as the cutoff distance (Bandyopadhyay, Valavala, & Odegard, 2009). Li and Strachan studied the process of polymerization and crosslinking in polymers through EPON-862 and DETDA using the DREIDING force field, with Buckingham van der Waals interactions, as presented in Figure 35 (C. Li & Strachan, 2010). At the macroscopic scale, the mass density, glass transition temperature, thermal expansion coefficient and elastic modulus calculated numerically were in good agreement with experiments. Furthermore, Yang (2010) studied tensile deformation and failure mechanisms in epoxy/copper bimaterials by means of largescale MD simulation, using PCFF force field equations, as shown in Figure 36.
59

Figure 35. Epoxy before and after Crosslinking between EPON-8672 and DETDA (Li and Strachan, 2010)
Figure 36. Snapshots of atomic configurations during pullout tests simulated with MD for copper-epoxy interfaces (Shaorui Yang, 2010)
60

3.1.2 MD model of an interface of silica and highly cross-linked epoxy resin (DGEBA and DETA)
We start our molecular investigation of epoxy by modeling the interface between silica, the main chemical component of Portland cement and granite, and a common epoxy resin used for structural components and coatings, Bisphenol A diglycidyl ether (DGEBA), mixed with a common curing agent, diethylenetri7mine (DETA). The molecular structures of DGEBA and DETA are illustrated in Figure 37.
Figure 37. Molecular structures of DGEBA (top) and DETA (bottom). Grey atoms are carbon, white atoms are hydrogen, red atoms are oxygen, blue atoms are nitrogen.
3.1.2.1 Simulation of cross-linking and equilibration We use the MD software MAPS to generate the initial DGEBA and DETA models,
with the stoichiometrically ideal composition ratio (5:2). The CH2 molecules are activated to react with DETA nitrogen atoms. Cross-linking, i.e. the process by with the hardener links with several epoxy molecules together towards a cured state, is simulated by
61

performing a series of multi-step relaxations that effectively reduce fluctuations in the creation of new bonds, which ensure a higher conversion rate (C. Li & Strachan, 2010).
All epoxy models are made of an amorphous cell containing a low density of resin and curing agent (0.5g/cm3), subjected to three-dimensional periodic boundary conditions. In order to study the sensitivity of our model to scale effects, we created two models of crosslinked epoxy, of two different sizes, as described in Table 6. By default, we use Model 2 for all the other sensitivity analyses and for the model calibration.

Table 6. Dimensions of cross-linked epoxy. W, D, H: model width, depth and height.

Model No.

Number of atoms

W x D ()

H ()

1

14,250

42 x 42

84

2

28,500

55.42 x 55.42

110.84

We use the DREIDING force field equation to account for inter-atomic forces within the epoxy bulk because of its excellent agreement with experiments reported in previous studies (C. Li & Strachan, 2010; Shaorui Yang, Gao, & Qu, 2013). The generated DGEBA and DETA box is equilibrated using the NVT ensemble simulation for 50ps (50 x 10-12 seconds) and then the NPT ensemble simulation for 400ps at atmospheric pressure, to ensure that the liquid reaches its equilibrium density. After this non-chemical equilibration step, simulations of chemical reactions are conducted. A chemical reaction occurs if the distance between two potentially reacting atoms is below a certain threshold AND if the reaction has priority in front of other potential reaction (e.g. if this is a primary reaction as opposed to a secondary reaction). The initial threshold distance used in our simulation is 5.64 (i.e., 5.64 x 10-10 m). A series of NPT simulations is performed to relax and
62

equilibrate the residual stress induced by the introduction of new chemical bonds in the previous steps. Figure 38 shows an epoxy molecule after cross-linking and Figure 39 shows an example of the cross-linked structure at 78% conversion rate (i.e. 78% cured, i.e. 78% cross-linked).
Figure 38. Epoxy after crosslinking between DGEBA and DETA (Arab, Shokuhfar, & Ebrahimi-Nejad, 2012)
63

Figure 39. MD specimen of highly cross-linked epoxy. D = W = 55 ; H = 110 .
3.1.2.2 Calibration of the MD epoxy model After cross-linking and equilibration, a series of simulations are conducted on cross-
linked epoxy to calibrate the epoxy glass transition temperature, coefficient of thermal expansion, mass density, elastic modulus, and Poisson's ratio of crossed-linked epoxy. The tensile strength and fracture energy are calibrated at the scale of the silica/epoxy interface, as explained in the next subsection.
Annealing simulations are first conducted to achieve the desired epoxy glass transition temperature and measure the coefficient of thermal expansion. An NVT ensemble simulation is conducted at 600K. Next, an annealing process is simulated to cool down the epoxy at room temperature. The NVT annealing process is simulated from 600K to 300K,
64

at the cooling rate of 10K/10ps. Nose-Hoover thermostat and barostat conditions are used. We use the DREDLING Buckingham (X6) form for van der Waals interactions. Results obtained during the cooling are plotted in Figure 40.

a

b

Figure 40. Properties of 78% cross-linked epoxy during the annealing simulation. (a)

Density vs. Temperature (b) Volume change vs. Temperature

A segmental linear regression method is utilized in order to determine the glass transition

temperature of the system. The predicted glass transition temperature of epoxy is in the

range 395-405K, which is similar to the glass transition temperature measured in

experiments reported in (Arab, Shokuhfar, & Ebrahimi-Nejad, 2012). The coefficient of

thermal expansion is calculated according to the following equation:

1 = W Y
( Q

where

L0

is

the

specimen

size

at

300K

and

OR OS

is the slope of linear of the box size versus

temperature. The coefficient is found equal to 714.28ppm/K in the rubbery state and

65

250ppm/K in the glassy state, respectively. This result is in the expected range of values (Shaorui Yang, 2010).
A NPT ensemble simulation is then run at 300K and atmospheric pressure in order to achieve the desired mass density of cross-linked epoxy at 300K. The calculated density at 300K is about 1.15 g/cm3, which is in good agreement with the experimental results presented in a previous study (1.16 to 1.19 g/cm3 according to Arab et al., 2012).
In order to calculate the elastic modulus and the Poisson's ratio of the cross-linked epoxy, we simulated uniaxial tension tests under the atmospheric pressure and a 300K temperature. The subsequent uniaxial tensile tests are conducted by using the MD software LAMMPS, which allows changing loading speed, conversion rate and geometry in the simulations. The cross-linked epoxy model constructed with MAPS is first imported into LAMMPS. A uniaxial strain is uniformly applied on the top of the cross-linked epoxy box. The bottom of the box is fixed to a layer of epoxy molecules. Figure 41 shows the simulated stress versus strain relationship of the 90% cured epoxy subjected to a uniaxial tension at a strain rate of 5 x 108 s-1. The elastic modulus and the Poisson's ratio of the 90% crosslinked epoxy are found to be in the range of 2.2-2.6 GPa and 0.35-0.3, respectively. These values are in agreement with those found in the experiments reported in (C. Li & Strachan, 2010; Xiong, Liu, Li, & Xie, 2002). Note that these experiments were conducted in quasistatic conditions, in contrast with our simulations, which were conducted at loading rates typical of MD simulations.
66

a

b

Figure 41. Simulation of a uniaxial tension test on 90% cross-linked epoxy at 300K,

atmospheric pressure and using a strain rate of 5 x 108 s-1. (a) Stress-strain response (b)

Lateral vs. longitudinal strain curve.

Sensitivity analyses were finally conducted to check if there is any dependence of the MD model of cross-linked epoxy to parameters other than the ones calibrated. We first study the influence of the loading strain rate. Axial stress-strain curves obtained at three loading rates are shown in Figure 42. Although the simulations indicate that the yield point, i.e. the point where the non-linear behavior initiates, strongly depends on the loading rate, the peak stress and the stress-strain curves before the yielding point do not depend on the loading rate in our simulations. We thus confirm that slow-rate experiments can be reproduced by high-rate MD simulations, with a satisfactory prediction of the mechanical response of the cross-linked epoxy. This is a significant advantage, because the computation time necessary to simulate epoxy tension tests in real time would make the study infeasible. In the following, all the simulations are conducted with a loading strain rate of 5 x 108 s-1.
67

Figure 42. Stress-strain response of the 90% cross-linked epoxy subjected to a uniaxial tension test at 300K, under different strain rates
The influence of the conversion ratio is illustrated in Figure 43, in which the silica/epoxy specimen is poxy specimen is subjected to a uniaxial tension test at a rate of 5 x 108 s-1 under a temperature of 300K, for 50%, 78%, and 90% cross-linking. As expected, the yield strain and the peak stress increase with the conversion rate. For instance, for tensile strains of 20% or more, the stress in 90% cross-linked epoxy is twice that of the 50% cross-linked epoxy and 1.5 times that of the 78% cross-linked epoxy.
68

Figure 43. Stress-strain response of cross-linked epoxy under a uniaxial tension test imposed at a rate of 5 x 108 s-1 under 300K, for different conversion ratios
Size effects are finally investigated by comparing the stress-strain curves obtained for the uniaxial tension test conducted at a rate of 5 x 108 s-1 under 300K on 90% crosslinked epoxy. Figure 44 shows the results obtained with models 1 and 2 described in Table 6. For tensile strains of 20% or more, the stress in model 1 (small size) is 80% of that in model 2 (medium size). The three sensitivity analyses indicate that the yield strain and peak stress are mostly influenced by the conversion ratio. The sample size has a lighter influence on the results and the strain rate exhibited almost no influence on the results.
69

Figure 44. Stress-strain response of cross-linked epoxy under a uniaxial tension test imposed at a rate of 5 x 108 s-1 under 300K, for two different sample sizes
3.1.2.3 Simulation of the silica/epoxy pull out test We now evaluate the mechanical behavior of the interface between epoxy and
concrete at the molecular scale by simulating a pullout test. Concrete is represented by a layer of crystalline silica. The unit cell of silica oxide was cleaved in the [0, 0, 1] direction on the interface side. The supercell of SiO2 is first generated in an amorphous cell in MAPS. The cell is then packed with DGEBA and DETA. Following the method presented in the previous subsection, cross-linking is simulated so as to achieve a conversion rate of 90%. The bi-material interface between silica and crosslinked epoxy is shown in Figure 45. Table 7 provides the dimensions of the interface
70

MD model. The layer of silica at the bottom is fixed. The top of the epoxy layer is subjected to a tensile strain, at a rate of 5x108s-1, to simulate the silica/epoxy debonding.

H

D

W

Figure 45. Geometry of the interface between highly cross-linked epoxy and silica

Table 7. Dimensions of the specimen used to simulate the pullout test at the interface between silica and highly cross-linked epoxy. W, D, H: width, depth, height.

Number of atoms

W x D ()

H ()

15,705 35.73 x 40. 3276

114.4

Pullout test results are shown in Figure 46. The stress-strain curve exhibits multiple peaks followed by stress relaxation every time debonding occurs. Figure 46 shows the influence of the conversion ratio on the interface response. The peak adhesive stresses at
71

the interfaces with 90% cross-linked epoxy are in average twice larger than those at the interfaces with lower epoxy conversion ratios.
Figure 46. MD simulation results obtained after simulating a pullout test on an interface between silica and 90% cross-linked epoxy. Top: Stress-strain curve. Bottom: Atomic configurations at key points of the stress-strain curve.
72

Figure 47. Stress-strain curves obtained after simulating a pullout test on an interface between silica and cross-linked epoxy, for different conversion ratios.
Our MD results are consistent with experimental and numerical data on DGEBADETA epoxy published in the literature, and show that MD can be used to explore the debonding properties of interfaces between concrete and other types of epoxy. In the following, we study the sealant used by G-DOT: High-Molecular-Weight-Methacrylate (HMWM).
3.1.3 Molecular models for the components of the concrete/HMWM interface Now that we demonstrated the feasibility of modeling concrete/epoxy interfaces
with MD, we now present an MD model to understand the fundamental processes that govern the mechanical behavior of interfaces between HMWM and concrete minerals (calcite and silica) at atomistic scale. We first present our MD models of calcite, silica and HMWM. Then, we present the MD models of mineral/HMWM interfaces and we describe the results of pullout tests that are simulated under different strain rates, temperatures and moisture conditions. We discuss the effect of strain rate last.
73

We first simulated the transition of HMWM to the glassy state, i.e. the process by which the hardener links with methyl methacrylate (resin) towards a cured state, also known as cross-linking. Figure 48 shows the molecules of resin and hardener, and Figure 49 illustrates the principle of the crosslinking reaction. We assumed that there was a freeradical polymerization of Methyl methacrylate (MMA). The radical initiator (cumene hydroperoxide) was decomposed to generate free radicals and to react with the MMA monomer (Yamamoto et al., 2013). The stoichiometric mixing ratio of methyl methacrylate and cumene hydroperoxide was 35:1, with a low initial density of 0.5 g/cm3. We used a 3D periodic boundary condition. The positions of the atoms were iteratively adjusted so as to minimize the potential energy of the HMWM system. Then, following previous studies (Hossain et al., 2010), the cross-linked system was subjected to an NVT ensemble simulation (imposed number of particles, imposed volume and imposed temperature) under 600K, for 5 x 10-9 s (105 timesteps of = 0.5). The HMWM crosslinked structure obtained by MD is a relaxed polymer network (traction free condition), as shown in Figure 50.

(a)

(b)

Figure 48. Molecular structure of (a) Cumene hydroperoxide (CHP), (b) Methyl methacrylate (resin)
74

3.1.3.1 Molecular models of HMWM and concrete minerals N
Polymerization

Methyl methacrylate

Polymethyl methacrylate

R
Methyl methacrylate Cumene hydroperoxide Figure 49. Schematic diagram of the chemical reactions involved in HMWM crosslinking.

75

Figure 50. HMWM atomic configuration after crosslinking. Carbon atoms are grey, oxygen atoms are red, hydrogen atoms are white.
Concrete mechanical properties are highly dependent on concrete structure and chemical composition, as highlighted in previous MD studies (Gao, 2018 and Xu, 2016). The most abundant material components in concrete and cement are quartzite and limestone. Limestone aggregates mostly consist of calcium carbonate (calcite), and sand and granite are mainly made of quartzite (silica). Therefore, we modeled calcite and silica at the molecular scale to investigate the resistance of concrete/HMWM interfaces to debonding. Figure 51 shows the lattices of calcite and silica built in MD. For calcite, the dimensions of the lattice are a = b = 4.990 , c = 17.061 , = = 90 = 120. For silica, the dimensions of the lattice are a = b = 4.913 , c = 5.405 , = = 90 = 120.
76

$

!
$

%

&
%

"

#

"

! &
#

Figure 51. MD lattice models of representative concrete minerals. For calcite, calcium atoms are orange, carbon atoms are black and oxygen atoms are red. For quartz, silicon atoms are red, oxygen atoms are white. Calcite: a = b = 4.990 , c = 17.061 . Silica: a
= b = 4.913 , c= 5.405

3.1.3.2 Material properties of HMWM After cross-linking and equilibration, we conducted a series of NPT simulations
(imposed number of particles, pressure and temperature) to calibrate the following properties: epoxy glass transition temperature, mass density, Young's modulus, Poisson's ratio. To perform these simulations, we used the DREIDING Buckingham (X6) form for the van der Waals interactions, and Nose-Hoover temperature thermostat and pressure barostat. In order to achieve a relaxed amorphous polymer structure, a temperature of 600K was first maintained at atmospheric pressure for 2.5 x 10-10 s (5 x 105 timesteps of = 0.5 ). Then, we simulated annealing with the NPT ensemble, by decreasing the temperature from 600K to 300K, within 2.5 x 10-10 s (5 x 105 timesteps of = 0.5).
77

Results obtained during the cooling stage are plotted in Figure 52. We employed a piecewise bilinear regression method in order to determine the glass transition temperature of HMWM, i.e. the temperature at the transition between the glassy and rubbery states of HMWM. The glass transition temperature predicted by MD is around 407K, which is conform to the glass transition temperature measured in published experiments (D. Meggers, 1998). Lastly, we simulated a relaxation test: a room temperature of 300K was applied for 2.5 x 10-10 s (5 x 105 timesteps of = 0.5). The calculated mass density at 300K was 1.053 0.5 g/cm3, which is in good agreement with the experimental results presented in a previous study (0.995 to 1.068 0.5 g/cm3) conducted by D. Meggers, 1998.
78

Volumetric change

Density [g/cm3]

0.25 0.2
0.15 0.1
0.05
!" = $%&'
0 250 300 350 400 450 500 550 600
Temperature [K]
1.1
!. #$%&/()% at 300K
1.05 1
0.95 0.9
0.85 250 300 350 400 450 500 550 600
Temperature [K]
Figure 52. Simulation of an annealing test on cross-linked HMWM, at atmospheric pressure, decreasing the temperature from 600K to 300K. Top: Volumetric deformation
vs. Temperature. Bottom: Density vs. Temperature
79

To calculate the tensile modulus and the Poisson's ratio of the cross-linked epoxy, we simulated uniaxial tension tests under the atmospheric pressure, at a temperature of 300K. The cross-linked epoxy model constructed with MAPS was first imported into LAMMPS. A uniaxial strain was uniformly applied at the top and at the bottom of the cross-linked HMWM model (called "box"). We first simulated uniaxial tension at a strain rate of 5 x 108 s-1, following previous work (Li and Strachan, 2010). Figure 53 shows the simulated stress versus strain relationship. For strains in the range 5-8 %, the tensile modulus and the Poisson's ratio of the cross-linked HMWM are found to be in the range of 2.076-2.7 GPa and 0.3-0.35, respectively. Published experimental data on HMWM mechanical properties were obtained by subjecting specimens to quasi-static loading under ambient temperature (D. Meggers, 1998). It is not feasible to model these experiments with MD, because of the computational cost of low strain-rate simulations with non-periodic boundary conditions. For example, for a strain rate of 1 x 108 s-1, a single uniaxial tension simulation took about 40 hours on Georgia Tech's Partnership for an Advanced Computing Environment (PACE), using 64 cores and 64 GB RAM. In MD, the response of large systems is typically simulated by applying periodic boundary conditions, which makes it challenging to capture strain localization and brittle failure. For instance, it is well-known that for polymer systems well below the glass transition temperature, MD simulations overpredict strength and do not capture brittle failure (Li and Strachan, 2010). It was also noted that compression strength of polymers increases with the loading rate (Richeton and Ahzi, 2007).
80

Stress [MPa]

200 Longitudinal stress Lateral stress
150

100

50

0

-50

-5

0

5

10

15

20

25

Strain [%]

10

Lateral strain [%]

8

6

4

2

0

0

5

10

15

20

25

30

Longitudinal strain [%]

Figure 53. Simulation of a uniaxial tension test on cross-linked HMWM at 300$K$, at atmospheric pressure and using a strain rate of 5 x 108 s-1. Top: Stress-strain response.
Bottom: Lateral vs. longitudinal strain curve.

81

We simulated the uniaxial tension test at three different strain rates (1 x 108 s-1; 5 x 108 s-1 ; 1 x 109 s-1). Figure 54 shows the stress-strain curves that we obtained numerically. We verify that the stress at the yield point increases with the loading strain rate (here, we define yield point as the point where the non-linear behavior initiates). We also note that, in average, when strains are in the range 4 - 7 %, the secant stiffness (Young's modulus) predicted numerically increases with the strain rate. This trend is consistent with the results obtained with MD by other authors (e.g., Richeton and Ahzi, 2007: Li and Strachan, 2010). In the discussion that follows, we use Richeton's model (Richeton and Ahzi, 2007) to explain how numerical predictions can be extrapolated to predict the HMWM tensile modulus at lower strain rates. In the following, all the MD simulations are conducted with a loading strain rate of 5 x 108 s-1, unless otherwise stated.
200

Stress [MPa]

150

100

50

Strain rate-1e9

Strain rate-5e8

Strain rate-1e8

0

0

5

10

15

20

25

Strain [%]

Figure 54. Effect of strain rate on the mechanical behavior of HMWM.

82

3.1.4 Molecular model of the mineral-HMWM interface This section aims to understand the fundamental processes that govern the
mechanical behavior of interfaces between HMWM and concrete minerals (calcite and silica) at the atomic scale, from an MD perspective.

3.1.4.1 Construction of the MD interface models We simulated pullout tests on HMWM/concrete interfaces at the molecular scale.
Concrete was represented by a layer of either calcite or crystalline silica. The lattices of calcite and silica are shown in Figure 55. Both mineral structures were cleaved in the [0, 0, 1] direction (on the interface side) to build the mineral substrate. The cleaved mineral blocks were generated in an amorphous cell in MAPS. A vacuum cell was added on top. The vacuum cells were then packed with crosslinked HMWM (Figure 49). The dimensions of the simulation's cells (W x D x H) were chosen based on previous studies (Xu and Wang, 2016; Yangming Gao et al, 2018), see Table 8.

Table 8. Dimensions of the simulation cells. W, D, H: width, depth, height.

Minerals

Number of atoms

W x D x H ()

Calcite

9,232

39.92 x 34.12 x 69.14

Silica

9,319

35.58 x 38.46 x 70.44

Geometry optimization and NPT equilibration tests were conducted with the interface models, in a similar way as for the HMWM model presented in the previous section. The interface debonding tests were simulated by fixing the bottom layer of the mineral substrate and by applying a tensile strain rate at the top of the HMWM layer. The

83

DREIDING force field was applied and the non-bonding forces were calculated with a cutoff distance of 12.

!

!

"

#

"

#

(a)

(b)

Figure 55. Molecular interface systems for (a) the calcite-HMWM model (b) the silicaHMWM model.

3.1.4.2 Results of the numerical interface debonding tests with different mineral substrates We first compared the results obtained with calcite and with silica for a debonding
test conducted at a rate of 5 x 108 s-1 under a temperature of 300K. A single pullout simulation took about 65 hours on Georgia Tech's Partnership for an Advanced Computing Environment (PACE), using 64 cores and 64 GB RAM. The fluctuations of the simulation results were smoothed using the moving average function implemented in MATLAB. The resulting stress-strain curves are shown in Figure 56. Both calcite and silica interfaces initially present a linear stress-strain response up to the peak stress, which is followed by a
84

significant stress drop, indicating a reduction of HMWM/minerals interactions. After debonding of the HMWM/mineral interface, stress converges to zero. The interfacial strength is higher for the HMWM/silica interface than for the HMWM/calcite interface.

200 Calcite Silica
150

Stress [MPa]

100

50

0

0

10

20

30

40

50

Strain [%]

Figure 56. HMWM/mineral interface debonding tests simulated at 300K, at a rate of 5 x 108 s-1, for two different concrete minerals silica and calcite.

We repeated the debonding tests in the presence of moisture. We added 100 water molecules in the interface system molecules by following the method by Y. Gao and collaborators (Yangming Gao et al, 2018). Starting with the interface model without water molecules, a thin confined layer of 100 water molecules was inserted between the mineral substrate and HMWM. The interface models with water molecules were then subjected to geometry optimization and equilibration procedures in a same manner as for interface models without water molecules. Stress-strain curves and snapshots of atomic configurations are shown in Figures 57 and 58. Stretching (phase A) is followed by microvoid formation (sequence B-C). Final separation is captured (phase D) and the critical

85

failure plane is found to be at the interface between HMWM and the mineral for both calcite and silica.
100

80

Stress [MPa]

60

40

20

$

!

"~#

0

0

10

20

30

40

50

Strain [%]

!

"

#

$

Figure 57. Calcite-water-HMWM interface system during the interface debonding tests at 300K, at a rate of 5 x 108 s-1, with 100 water molecules.

86

Stress [MPa]

120

100

80

60

40

20

!

"~#

$

0

0

10

20

30

40

50

Strain [%]

!

"

#

$

Figure 58. Silica-water-HMWM interface system during the interface debonding tests at 300K, at a rate of 5 x 108 s-1, with 100 water molecules.

87

3.1.4.3 Results of the numerical interface debonding tests under different temperature and

moisture conditions

We studied the response of the HMWM/water/mineral interface systems at four temperatures (200 K, 250 K, 300 K and 400 K). The strain rate was 5 x 108 s-1 and 100

water molecules were present at the interfaces. The simulation results indicate that the

primary separation occurs at the interface between HMWM and the mineral substrate. The

deformation of the HMWM bulk does not exhibit any dependence to temperature. Figure

59 shows that the strength of both the HMWM/silica and HMWM/calcite interfaces

decreases when temperature increases. For example, for the HMWM/calcite interface

system, we note a 16% decrease of interface strength from 200 K to 300K and a 63%

decrease of interface strength between 300 K and 400 K. The dependence of tensile

strength to temperature was expected for this range of temperature, since the glass

transition temperature is around 407 K.

200

calcite

150

silica

Strength [MPa]

100

50

0 200 250 300 350 400
Temperature [K]

Figure 59. Interface strength predicted by MD under different temperatures.

88

Interface strength and work of separation were calculated for both substrates during

the interface debonding tests performed at 300 K, with and without moisture. Results are presented in Figure 60. The work of separation was calculated as the area of the stressdisplacement curve divided by the energy required to separate a unit area of interface. The

strength of the silica/HMWM interface is 30 MPa higher than the HMWM/calcite interface, regardless of moisture conditions. In the presence of moisture, interface strength decreases by 32% for the calcite/HMWM interface and by 26% for the silica/HMWM interface. In dry conditions, the work of separation is about 175mJ/m2 for both the calcite and the silica systems. In the presence of moisture, the work of separation is 125.97mJ/m2 for the calcite/HMWM interface and 156.69mJ/m2 for the silica/HMWM interface.

Dry condition 200

Wet condition

dry condition 150

Wet condition

Work of separation [mJ/m2 ] Strength [MPa]

160

120

120

90

80

60

40

30

0
calcite

quartz

0
calcite

quartz

(a) Tensile strength

(b) Work of separation

Figure 60. Interface strength predicted by MD under different moisture conditions.

89

3.1.4.4 Results of the numerical interface debonding tests under different strain rates We simulated the interface debonding test at three strain rates: 1 x 108 s-1; 5 x 108 s-1
; 1 x 109 s-1, under a 300 K temperature. The stress-strain curves obtained for the calcite and silica interface models are shown in Figure 61. The shape of the stress-strain curves and the peak strength are influenced by the strain rate. We note that the peak stress is higher than the peak stress reported in experimental studies (D. Meggers, 1998). This is because the strain rates used in the MD simulations were way higher than the rates used in the experiments. As explained above, the computational cost of the MD simulations is prohibitive, which does not allow analyzing the sensitivity of our interface model at quasistatic strain rates. We will discuss how to overcome this limitation in the following.

Stress [MPa] Stress [MPa]

150

200

Strain rate-1e8 Strain rate-5e8

Strain rate-1e8 Strain rate-5e8

Strain rate-1e9

150

Strain rate-1e9

100

100

50 50

0

0

10

20

30

40

50

Strain [%]

0

0

10

20

30

40

50

Strain [%]

(a) Calcite substrate

(b) Silica substrate

Figure 61. Effect of strain rate on the mechanical behavior of the HMWM/mineral interface during the interface debonding tests.

90

The simulation results indicate that for both mineral substrates, the higher the strain rate, the higher the peak strength and the higher the total area of the stress/strain curve. These results are consistent with observations made by previous researchers (Shaorui Yang et al., 2013). We also note that the elastic behavior of the interface systems does not exhibit significant dependency to strain rate.

3.1.4.5 Interface strength and fracture energy

The stress-strain response predicted in the MD interface tension tests is similar to a

cohesive law, which led some authors to use the MD stress-strain curve as a traction-

separation law to define a Cohesive Zone Model (CZM) at the metric or submetric scale

(A. Sazgar et al., 2017; Xu and Wang, 2016; Wu Gui Jiang et al., 2018). Figure 62

illustrates this method. As an example, we fitted an exponential traction-separation law (in

red) to the stress-strain response obtained during the interface debonding test conducted

with the dry HMWM/silica interface system at a loading rate of 1 x 108s-1, under a

temperature of 300 K and at atmospheric pressure. The exponential softening curve has the

following equation:







=

6

WY 6



W1

-

Y 6

where 6 is the peak stress and 6 is the critical strain, defined as the strain that is reached

when = 6.

91

1

Figure 62. Simulation result and curve fitted exponential traction-separation law (in red) to the stress-strain response obtained during the interface debonding test conducted with
the dry HMWM/silica interface system.
After fitting, it is found that the strength of the dry interface is 152.22MPa and the critical strain is 9%, with an R-square error of 95%. As shown in Figures 57 and 58, the interfacial strength and the critical strain of the HMWM/water/calcite system are 66.81MPa and 7.5% respectively, while the interfacial strength and the critical strain of the HMWM/water/silica system are 86.89MPa and 7.34% respectively. These interfacial strengths are 16 times larger than the interfacial strengths found experimentally (Transpo Inc., 2017). This difference is due to size effects and to strain rate effects. Size effects cannot be avoided in that sense, we conclude that the stress-strain curve obtained by MD cannot be used as is in a macroscopic CZM model. However, the Richeton's model can be calibrated to extrapolate the MD predictions made at high strain rate to lower strain rates. We calibrate the Richeton's model in the following section.
92

3.1.4.6 MD predictions at low strain rate with phenomenological models

The effects of strain-rate and temperature on the mechanical behavior of composites

and polymers were investigated by using the Richeton's model and the Johnson-Cook

model's (Shahin Shadlou et al., 2014; T J Jolmquist et al., 2016; Zhemin Jia et al., 2016;

Huiping wu et al., 2019). In this study, we used these two phenomenological models to

predict the tensile modulus of HMWM and the strength of the HMWM/concrete interface.

Based on the Richeton's model, the elastic modulus of polymers is expressed as (Richeton

and Ahzi, 2007):

T

=

T ;>9

L1

+



W;> 9YM

where T and T;>9 are the elastic modulus and the reference elastic modulus under a

static loading rate, S is strain rate sensitivity constant, and ;> 9 are the effective plastic

strain rate and the reference strain rate. We calibrated the S parameter to fit the Richeton's

model to the MD results of HMWM uniaxial tension tests. For the reference elastic

modulus, the tensile modulus of HMWM is 0.28GPa at a tensile strain rate of 0.1s-1,

according to experimental work reported in D. Meggers,1998. We found a good fit between

the Richeton's model and the MD results with S=0.3428. The predictions of the Richeton's

model, shown with a red dashed line in Figure 63, could be used to predict the values of

the tensile modulus under low strain rate if more experimental results were available to

check the extrapolation.

93

Tensile modulus [GPa]

3.5

Richeton model

3

Experiment

MD simulation

2.5

2

1.5

1

0.5

0 100

105

1010

Strain rate [/s]

Figure 63. Variation of the tensile modulus of HMWM as a function of strain rate in a uniaxial tension test. Experimental data is used as reference. Blue dot: reference experimental data. Black dots: MD results.

The Johnson-Cook model expresses the dependence of interface tensile strength to

the strain rate and to temperature, according to the following equation (Johnson and Cook,

1985}:

T

=

T ;>9

1

+



W;> 9Y

J1

-

L



-

;<<+

+
M

N

U - ;<<+

where T and T;>9 are the yield strength and the reference yield strength under a static

loading rate, C is a strain rate sensitivity constant, ;<<+ and U are the room temperature

and the glass transition temperature, m is a temperature sensitivity parameter. Based on the

PCC-SSD bond strength test results reported in the technical datasheet of Transpo company

(Transpo Inc., 2017), the tensile strength of the HMWM/concrete interface is 4.2MPa at a reference strain rate of 0.1 s-1. We calibrated the values of parameters C and m to match

94

the MD results of interface debonding tests. We found C=1.019 and m=0.441. The predictions of the Johnson-Cook model, represented by a hyper-surface in Figure 64, could be used to predict the variations of interface tensile strength with temperature and strain rate, if more experimental data was available to check the extrapolation. In general, interfacial strength increases when temperature decreases and when the strain rate increases.

MD simulation Experiment Johnson-cook model

Interfacial strength [MPa]

150

100

50

0 200
250 300 350
Temperature [K]

400 1010

100 105
Strain rate [/s]

0

20

40

60

80

100

120

140

Figure 64. Interfacial strength as a function of strain rate and temperature using Johnson-Cook model for the HMWM/calcite interface model. The red dot is an experimental data point. The blue dots are the data points obtained by MD simulation.

95

3.1.4.7 Contribution of the non-bonded energy to the work of separation of the

HMWM/mineral interfaces

Figure 65 presents the evolution of the work of separation, the non-bonding energy

and the bonding energy of the HMWM/mineral interfaces during the interface debonding

tests simulated by MD. The non-bonding energy Enb and the bonding energy Eb were

calculated according to DREIDING force field equation. Results show that the work of

separation is always almost equal to Enb+Eb. We note that the bonding energy oscillates

between positive and negative values. This trend was already noted in (Xiaoying Zhuang,

2018).

400

300

Total potential energy Non-bond energy

Total potential energy Non-bond energy

300

Bond energy 200

Bond energy

Energy [kcal/mol] Energy [kcal/mol]

200

100

100

0

0

-100

0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

Strain [%]

Strain [%]

(a) Calcite substrate

(b) Silica substrate

Figure 65. Energy distribution during the interface debonding tests at 300 K for a strain rate of 5 x 108s-1 under the dry condition.

Additionally, the contribution of the bonding energy to the work of separation is negligible for the HMWM/silica interfaces. The bonding energy was equal to a fifth of the nonbonding energy for the HMWM/calcite interfaces. For both HMWM/mineral interfaces,
96

van der Waals forces were thus the main forces that resisted debonding. This somewhat surprising result was already highlighted in previous studies, e.g. Shaorui Yang et al., 2013; Xu and Wang, 2016; Yangming Gao et al., 2018; Wu Gui Jiang et al., 2018; Xiaoying Zhuang, 2018.
3.1.5 Summary and conclusions We analyzed the mechanical properties of DGEBA-DETA epoxy/silica and
HMWM/concrete interfaces based on Molecular Dynamics (MD). The MD models of epoxy and HMWM are built in several phases, including geometry optimization, heating, annealing and relaxation. The glass transition temperature and the mass density predicted by the MD models are conform to experimental data reported in the literature for both DGEBA-DETA epoxy and HMWM. Uniaxial tension tests showed that all interfaces fail by debonding, at the surface of contact between epoxy (DGEBA-DETA or HMWM) and the mineral substrate, and that the interfacial strength decreases in the presence of moisture, under low strain rate, at high temperature, or at low conversion (crosslinking) ratio. MD results showed that the work of separation was mostly done by van der Waals forces on HMWM/concrete minerals. We used published experimental data at low strain rate along with our MD results at high strain rate to calibrate Richeton's model and Johnson-Cook model. Richeton's model relates the tensile modulus of HMWM to the strain rate. Its successful calibration indicates that our MD results could be extrapolated to predict the tensile modulus of HMWM at low strain rate. Johnson-Cook model predicts interfacial strength as a function of temperature and strain rate, given one reference interfacial strength at given temperature and strain rate conditions. MD results confirm that HMWM should be applied on dry surfaces and in concrete exposed to mild temperatures. Additionally, MD
97

results suggest that HMWM application in concrete with high silica contents is more likely to last than in concrete with high calcite contents.
3.2 Simulation of concrete/HMWM at metric scale This section provides a quantitative assessment of concrete repaired by injection of
High Molecular Weight Methyl methacrylate (HMWM). We simulated the experiments with the Finite Element Method which are subjected to uniaxial compression and splitting tension tests with and without reparation. 3.2.1 HMWM model for FE simulation
To simulate the tests conducted in Phase III on repaired concrete, we used the CDP model for non-damaged concrete, and we modeled the cracks filled with cured epoxy (HMWM in the glassy state) with bilinear softening cohesive zone elements (CZE) that capture both tension and shear strength, as illustrated in Figure 66.
!
#$
" Figure 66. HMWM linear softening traction-separation law
98

The traction-separation law used in ABAQUS is expressed in terms of nominal

tractions and displacement, as follows:

?

?? 0 0 ?

= ,,= ... = 0 == 0 ^ ,,= ... =

-

0 0 -- -

T is the traction force vector, with components in the normal direction (?) and in

the shear (tangential) directions =and -). The cohesion stiffness tensors are noted 88 ,

where i stands for n, s or t. is the displacement vector, corresponding to nominal strains.

We

assigned

??

a

value

50

times

(VP
-

=

50)

larger

than

that

of

the

Young's

modulus

(which was given as 0.95GPa by D. Meggers, 1998), and == and -- are obtained by shear

modulus formula proposed as (Jin and Arson, 2017):





?? = , == = -- =

In the simulations, we used the strength parameters provided by the supplier (11 MPa

in tension, 14 MPa in shear). We used a power law (J. Edward et al., 1965) to relate the

critical tensile and compressive fracture energies and account for mixed-mode CZE

debonding, as follows:

L??W

V
M

+

L==W

V
M

+

L--W

V
M

=

1

where ?W , =W , and -W are the critical fracture energies to initiate fracture

propagation in the tension, first shear, and second shear directions, respectively. The

power-law parameter was calibrated by simulating the uniaxial compression tests done

on repaired concrete cylindrical specimens, as explained in the following section. To ensure stability, we used a time increment smaller than the following critical time increment (Dassault systems, 2012):

99



=

6

6 S6

where 6 is the thickness of the cohesive layer, 6 is the elastic modulus of HMWM and 6 is the density of HMWM.

3.2.2 Verification of the CDP model for concrete repaired by epoxy The material parameters that could not be found in the literature or in the
specifications given by the suppliers were calibrated to ensure that the force-displacement curves obtained numerically fit the experimental results. Note that in the simulations of repaired reinforced concrete (RC) beams presented in Chapter 4, we used repaired concrete Finite Elements (FE) in addition to the CZEs. Table 9 summarizes all the material parameters used in the FE model. In the following, we explain our calibration simulations.

100

Table 9. Material parameters used in the FEM model of repaired concrete.

Plain and repaired concrete elements

Poisson's ratio

Standard value

6(

Literature review



Literature review

Dilation angle

Literature review

Eccentricity

Literature review

G(/6

Literature review

6

Literature review

Plain concrete elements

Young's modulus

From experiments

Tensile strength

From experiments and equation

Fracture energy

From experiments and equation

Repaired concrete elements

Young's modulus

From repaired compression tests

Tensile strength

Calibrated

Fracture energy

Calibrated

Cured HMWM Cohesive Zone Elements

Young's modulus

Literature review

Poisson's ratio

Standard value

Cohesive stiffness

Standard practice

Tensile strength

Calibrated

Shear strength

Literature review

Tensile fracture energy

Calibrated

Shear fracture energy

Calibrated

Power law

Calibrated

CZE thickness

Calibrated

Steel truss elements

Young's modulus

Standard value

Yield tensile strength

Standard value

Ultimate tensile strength

Standard value

Ultimate tensile strain

Standard value

0.2 0.02 0.2 30 0.1 1.12 0.75
28.52 GPa 2.28 MPa 0.105 N/mm
31.82 GPa 7.6 MPa 0.3 N/mm
0.95 GPa 0.3
47,500 GPa 8 MPa 14 Mpa
0.4 N/mm 1.2 N/mm
1 0.15 mm
200 GPa 420 MPa 520 MPa
0.1

101

3.2.2.1 Uniaxial compression tests on repaired concrete cylinders In Phase III, Uniaxial compression (UC) tests were conducted on cylinders 150mm
in diameter and 300mm in height. Concrete cylinders were cured and aged in a controlled fog room for 28 days. The HMWM-repaired concrete cylinders were fabricated by cutting plain concrete cylinders lengthwise on a tile cutter. Once cut, the faces of the two halves were ground flat in order to have a smooth surface for bonding. Construction silicon caulk was then applied to the edges of the sides of the two halves and a slim cardboard spacer was placed at the top and bottom to create a thin space (``crack''). The two halves were then pushed and tightened together using hose clamps, as illustrated in Figure 67. The silicon caulk was allowed to cure for 24 hours. The epoxy was poured into the top opening of the sample until the newly constructed ``crack'' was completely filled with HMWM. The repaired concrete specimens were cured for at least 48 hours before testing.
Figure 67. HMWM-repaired concrete specimens used in the uniaxial compression tests UC tests were conducted according to the ASTM C39 standards, for both the plain
and repaired cylinders. A SATEC PRISM MKIII-C with the INSTRON 59-R7 controller was used to apply loading at 60 kips (420 MPa) per minute. Three tests were conducted on HMWM-repaired concrete, and the average compressive strength was 40.54 MPa.
102

The repaired specimens were modeled in the same way as the plain concrete specimens studied in Chapter 2, except that HMWM CZ elements were inserted in the middle longitudinal plane. The rigid platens at the top and bottom of the specimen were modeled explicitly with rigid FEs. The interface between the specimen and the platen was modeled with a hard contact (non-penetrability condition), with a friction coefficient of 0.1. Fixed displacements were imposed to the centroid of the bottom platen. The vertical displacement at the centroid of the top platen was increased until the ultimate axial compressive displacement was reached. We extracted the force and the displacement at the centroid of the top rigid platen during the UC test simulations. Figure 68 (a) shows the vertical stress/displacement curves obtained numerically for the plain and repaired concrete specimens. For both types of specimen, the peak stress found numerically matches the tensile strength measured experimentally (32.54 MPa for plain concrete and 40.54 MPa for repaired concrete). Computational results indicate that the initial stiffness of plain concrete is lower than that of repaired concrete. Figure 68 (b) and (c) show a parametric study of the CZ thickness and of the exponent parameter . Based on these sensitivity analyses, we found that the best fit was 0.15 mm for the CZ thickness and 1 for the exponent . The value of is consistent with previous studies, which recommend a value of 1 for epoxy adhesive materials (Camanho and Davila, 2002). The distribution of damage in the repaired concrete specimen is shown in Figure 68 (d). The damaged zone obtained numerically is in good agreement with the fracture pattern observed experimentally during the UC tests: cracks appear mostly at the bottom of the specimen and along the reparation plane.
103

Stress[MPa]

Compressive strength [MPa]

60

50

Plain concrete

Experiment

50

Repaired concrete FE-Plain concrete

FE-Concrete cylinder

FE-Repaired concrete

45

40

30

40

20 35
10

0

30

0

0.5

1

1.5 0.05 0.1 0.15 0.2 0.25 0.3 0.35

Displacement [mm]

HMWM Thickness [mm]

(a)

(b)

Stress[MPa]

60 50 40 30 20 10

Experiment FE- :1 FE- :1.5 FE- :2

SDEG
(Avg: 75%)
+1.973e+00
+1.000e+00 +9.167e01 +8.333e01 +7.500e01 +6.667e01 +5.833e01 +5.000e01 +4.167e01 +3.333e01 +2.500e01 +1.667e01 +8.333e02
+0.000e+00 7.785e01

DAMAGEC
(Avg: 75%)
+2.022e+00
+1.000e+00 +9.167e01 +8.333e01 +7.500e01 +6.667e01 +5.833e01 +5.000e01 +4.167e01 +3.333e01 +2.500e01 +1.667e01 +8.333e02
+0.000e+00

0

0

0.5

1

1.5

Displacement [mm]

(c)

ODB: repairerate01.odb Abaqus/Explicit 6.121 Wed Nov 27 14:59:00 Eastern Standard Time 2019

Step: Step1

Y

X

Increment 157689: Step Time = 0.3000

Primary Var: SDEG

Deformed Var: U Deformation Scale Factor: +1.000e+00

Status Var: STATUS

Z

(d)

ODB: repairerate01.odb Abaqus/Explicit 6.121 Wed N

Step: Step1

Y

X

Increment 157689: Step Time = 0.3000

Primary Var: DAMAGEC

Deformed Var: U Deformation Scale Factor: +1.000e+00

Status Var: STATUS

Z

Figure 68. Comparison of numerical and experimental results for the UC test. (a) Stressdisplacement curves obtained numerically, and strengths obtained experimentally. (b) Numerical strength of the repaired specimen for different CZ thicknesses. (c) Numerical stress-displacement curves obtained with repaired concrete specimens, for different
values of the exponent . (d) Contour plot of damage at the end of the UC test simulation in repaired concrete, compared to a photograph of the fractures observed experimentally.

104

3.2.2.2 Notched repaired three-point bending tests without reinforcement In order to analyze the mode I interfacial debonding behavior of repaired concrete
and measure the fracture energy of concrete before and after reparation, we conducted three-point bending (TPB) tests, following Hillerborg's recommendations (A. Hillerborg, 1976), see Figure 69. The notch was created by inserting a single wall corrugated cardboard with a thickness of 3 mm and a depth of 82 mm. The HMWM-repaired beams were produced by cutting plain concrete beams in two halves with a masonry saw. A thin spacer was placed between the two halves, which were then reattached and sealed with silicon caulk. A corrugated cardboard 3 mm thick, 82 mm deep was placed in the thin space between the two halves. The silicon cured for 24 hours. Then HMWM was poured in the thin space and left in place to cure for 48 hours. Lastly, the corrugated cardboard was removed, to create a notch in the repaired beam.
Figure 69. Non-repaired (top) and repaired (bottom) non-reinforced concrete beams after three-point bending tests.
105

A 50 kips (345 MPa) Interface 1220AO-50k load cell was employed to apply a load at a rate of approximately 2 kips (14 MPa) per minute. The crack mouth opening displacement (CMOD) was measured by means of a lateral variable differential transmitter (LVDT). We used a RDP DCTH200AG LVDT, with a maximum travel distance of 5 mm. Two non-repaired beams and two repaired beams were tested. The peak force of the HMWM-repaired beams was higher than that of the non-repaired beams. In both the plain and repaired specimens, vertical cracks propagated at the bottom of the beam (see Figure 69). In the repaired beams, the cracks formed within the concrete, and not at the interface between concrete and HMWM.
We used a finer mesh size of 1 mm in the region of the expected crack path. The load/CMOD curves obtained numerically lie between the load/CMOD curves obtained experimentally, for both the plain and repaired beams (Figure 70). In other words, the model predicts a representative average behavior of plain and repaired concrete in tension. The behavior of the repaired beams was similar to that of the plain beams, with a slightly higher peak load (11.44 kN instead of 10.84 kN in the simulations). As expected, tensile damage propagates vertically ahead of the notch, within the concrete, along the crack plane filled with HMWM (i.e., along the CZ elements).
106

15

Load [kN]

10

5

Plain concrete

Repaired concrete

FE-Plain concrete

FE-Repaired concrete

0

0

0.05

0.1

0.15

0.2

CMOD [mm]

Figure 70. Evolution of the Crack Mouth Opening Displacement during the TPB tests done on plain and repaired non-reinforced concrete beams: numerical predictions are
representative of the average behavior observed experimentally.

3.2.3 Summary and conclusions Cut-and-sealed specimens were subjected to uniaxial compression tests and three-
point bending tests. The tests were repeatable and results indicated that repaired concrete specimens had higher stiffness and higher strength than plain concrete specimens. The experiments were simulated with the FEM. Concrete volume elements were assigned a constitutive model that combines continuum mechanics and plasticity, while the concrete/HMWM interface was modeled with CZ elements with a bilinear softening traction-separation law that allowed mode I, mode II, mode III and mixed mode interface debonding. Model parameters were calibrated against the uniaxial compression test results and the calibration was validated against the three-point bending tests.

107

CHAPTER 4. NUMERICAL SIMULATION OF REPAIRED RC AND PSC BEAMS
In this chapter, we present simulations of reinforced beams, with and without prestress, before and after reparation. We start by validating the CDP model enriched with CZ elements against results of three-point bending tests conducted on reinforced concrete (RC) beams. Next, we simulate stress and damage evolution during RC beam pre-stressing. We simulate a four-point bending test performed on a PSC beam by other authors and we compare our simulation results to those obtained if the pre-stress cracks are repaired before loading. Lastly, we analyze the sensitivity of PSC beams to initial crack distributions described in the PCI guidelines. We compare the responses of the pre-damaged beams to the responses of the same beams, after reparation. We compare the response of as-built and repaired PSC beams to that of RC beams of same volume fraction to distinguish the effects of rebar volume fraction vs. pre-stress on the performance of the reparation by HMWM deposition.
4.1 Three-point bending tests with RC beams In Phase III, the reinforced concrete (RC) beams were fabricated with the same
dimensions as the non-reinforced beams, without any initial notch. Two #6 steel rebars with a yield strength of 420 MPa were placed longitudinally. Two RC beams were tested as-built, to investigate the force-displacement behavior without HMWM injection. Three beams were loaded to develop different sizes of cracks (minor cracks obtained under a load of 80 kN or below, and major cracks obtained under a load of 100 kN). The three beams were then unloaded, removed from the load frame and laid on their sides. The change in orientation of the beam allowed gravity to assist the epoxy flow as it filled the cracks.
108

HMWM was cured for at least 48 hours before testing. The cured repaired beams were then tested to failure with the same method as for the TPB tests done on non-reinforced beams. The displacement fields measured by DIC were post-processed to calculate the strain fields. The post-processed field of maximum principal strain is shown in Figure 71 for beams with minor and major cracks.
4.1.1 Validation against experiments We modeled the RC beams with 3D solid concrete elements with an edge size of
10 mm and two-node linear truss steel elements with a length of 10 mm. We simulated a TPB test with the parameters given in Table 9. The load-displacement curve obtained numerically is compared to the load-displacement curves obtained experimentally (Figure 69.a). Each RC beam had a speckled longitudinal face that was filmed during the TPB test. We used Digital Image Correlation (DIC) to map the displacement field on those speckled faces. We then calculated the field of maximum principal strains. We found that cracks were present wherever the maximum principal strain exceeded 0.002. We used that threshold to define the damaged zone after the TPB tests, as illustrated in Figure 72. The contour of the damaged zone identified by DIC was then transformed into a binary image and meshed, by using the tracing image method in Solidworks. In that mesh, HMWM CZEs were inserted between all FEs, by using Truster's 3D interface element insertion code (T J truster, 2015). The volume elements were assigned the CDP model, with the properties of a "repaired" concrete, i.e. concrete with a larger stiffness, larger tensile strength and larger compressive strength than plain concrete (see Table 9). The "repaired concrete" elements represent concrete that contains small cracks filled with HMWM (i.e., cracks smaller than
109

the fractures represented by the CZEs). The model thus accounts for the propagation and repair of both large cracks (of the order of 0.01 m in width) and small cracks (less than 0.01 m in width).

(a)

(b)

(c)

(d)

Figure 71. Photographs of the RC beams during the TPB tests and the corresponding field of maximum principal strain obtained by DIC (a) when major cracks developed (loading up to 100 kN); (b) when the beam with repaired major cracks failed; (c) when minor cracks developed (loading up to 80 kN); (d) when the beam with repaired minor
cracks failed.

110

Figure 72. Method to identify the damaged zone after the TPB tests on RC beams. Top: the field of maximum principal strain is calculated by DIC. Bottom: the zone where the
maximum principal strain exceeds 0.002 is binarized and meshed with SolidWorks.

Out of the three repaired RC beams, two were loaded up to 100 kN and exhibited

"major cracks" before reparation, and one was loaded up to 80 kN and exhibited "minor

cracks" before reparation. We categorized cracks as minor or major depending on their

crack width, which were calculated based on the following equation (G. J. Schreppers et

al., 2011):

6;

=

GF

\X

-

X^

where E is the Young's modulus of concrete, 6; is the crack width and GF is the crack bandwidth, which is equal to the length of the edge of an element in the numerical model.

X and X are respectively the maximum principal strain and maximum principal stress,

which were post-processed from the simulation results. We found that the major cracks

(obtained when the beam was loaded up to 100 kN) had a width over 0.1 mm. For minor

cracks (obtained when the beam was loaded up to 80 kN), the main crack width was less
111

than 0.1 mm. We simulated the three-point bending (TPB) tests on the three beams after reparation. Figure 73 (b) shows that the numerical result predictions for repaired RC beams with major cracks are in good agreement with the experimental results up to the peak. After the peak, we got two very different experimental responses for the two beams with major cracks, and the numerical curve matches one of the two experimental curves. So overall, the numerical model performs satisfactorily. Figure 73 (c) shows that the proposed model reproduces well the behavior of the repaired RC beam with minor cracks, which is actually very similar to that of as-built RC beams. Simulations show that the capacity of the repaired RC beams with major cracks was approximately 40% higher than that of the as-built RC beams. This result suggests that epoxy injection could strengthen concrete with larger cracks with a width of 0.1mm and above, for which the specific area of contact between the concrete and HMWM is high.
The numerical model predicts crack or damage propagation in the zone where cracks were observed experimentally, and it predicts a maximum principal strain field in good agreement with that obtained by DIC (Figure 74). The model performs equally well for the specimens that had major cracks and for the specimens that had minor cracks. In particular, the model captures the occurrence of flexural and shear cracks at the same locations as in the experiments.
112

Load [kN] Load [kN]

150
100
50
0 0

250

RC concrete #1

RC concrete #2

FE simulation

200

Repaired after major crack #1 Repaired after major crack #2 FE simulation

150

100

50

2

4

6

8

Displacement [mm]

0

10

0

5

10

15

Displacement [mm]

(a)
150
100

(b)
RC concrete Repaired after minor crack FE simulation

Load [kN]

50

0

0

5

10

15

Displacement [mm]

(c)

Figure 73. Numerical and experimental load/displacement curves obtained during TPB

tests performed on (a) As-built RC beams; (b) Repaired RC beams with major cracks; (c)

Repaired RC beams with minor cracks.

113

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 74. Experiment and simulation results obtained for repaired RC beams: (a)

Photograph of repaired RC beam with major cracks; (b) Photograph of repaired RC

beam with minor cracks; (c) Snapshot of the FEM mesh, showing the zones of damage at

failure for the repaired beam with major cracks; (d) Snapshot of the FEM mesh, showing

the zones of damage at failure for the repaired beam with minor cracks; (e) DIC results

at failure, showing the field of maximum principal strain in the repaired beam with major

cracks; (f) DIC results at failure, showing the field of maximum principal strain in the

repaired beam with minor cracks; (g) Snapshot of the FEM mesh, showing the

distribution of maximum principal strain at failure in the repaired beam with major

cracks; (h) Snapshot of the FEM mesh, showing the distribution of maximum principal

strain at failure in the repaired beam with minor cracks.

4.1.2 Sensitivity analyzes We analyzed the sensitivity of the mechanical response of the repaired RC beams
to the width of the cracks filled by epoxy. Based on the FE simulation results on RC beams which were loaded at 80kN and 100kN, we first identified a target zone for repair, defined as the zone in which the crack width exceeded a certain threshold after the TPB test. We inserted CZ elements and replaced the concrete FEs by repaired concrete FEs in the target
114

zone, for different threshold values. For the RC beams loaded up to 100 kN, major cracks occurred. We compared the response of the RC beams where only cracks of width above 0.1 mm were repaired to RC beams where cracks of width above 0.05 mm or above 0.01 mm were repaired. Results of crack width above 0.01mm are presented in Figures 75 (a) and (b), where the grey zone in the mesh represents the target repair zone. We note that when cracks of width at least 0.01 mm are repaired, the model predicts the load displacement curve obtained experimentally. When the width threshold is set to a larger value, the model tends to under-predict the strength of the repaired beam.
For the RC beams loaded up to 80 kN, minor cracks occurred (maximum crack width: 0.036mm). When we simulated the TPB test with a width threshold of 0.1 mm, we verified that the response of the repaired RC beam was similar to that of the damaged asbuilt RC beam. This was expected, since the maximum crack width observed during the experiments at a load of 80 kN was 0.036mm < 0.1 mm. We then compared the response of the RC beams where only cracks of width above 0.01 mm were repaired to RC beams where cracks of width above 0.005 mm were repaired. Results are presented in Figures 75 (c) and (d). We note that when cracks of width at least 0.01 mm are repaired, the model predicts the load - displacement curve obtained experimentally.
To summarize, our simulation results suggest that HMWM can penetrate cracks of width 0.01 mm and above by gravity. Additionally, the results presented in the previous section indicate that when cracks of width 0.1 mm and above occur, repaired RC concrete has a strength that can exceed the strength of the as-built RC concrete by 30% to 40%.
115

Crack width 0.01mm (Cracks at 100kN)

Crack width 0.01mm (Cracks at 80kN)

(a)

300

Repaired after major crack #1

250

Repaired after major crack #2 Crack repaired > 0.01mm

Crack repaired > 0.05mm

200

Crack repaired > 0.1mm

(b)

300

RC concrete #1

250

RC concrete #2 Repaired after minor crack

Crack repaired > 0.01mm

200

Crack repaired > 0.005mm

Load [kN] Load [kN]

150

150

100

100

50

50

0

0

0

5

10

15

0

5

10

15

Displacement [mm]

Displacement [mm]

(c)

(d)

Figure 75. Simulation of the TPB test with RC beams, after repairing cracks of different

widths. (a) Maximum load 100kN, crack width threshold 0.01 mm; (b) Maximum load

80kN, crack width threshold 0.01 mm; (c) Load-displacement curve for RC beams loaded

up to 100 kN, with repaired cracks of width > 0.01 mm, 0.05 mm, 0.1mm; (d) Load-

displacement curve for RC beams loaded up to 80 kN, with repaired cracks of width >

0.005 mm, 0.01 mm.

116

4.2 Numerical simulation of repaired prestressed concrete girders In this section, we present a method to simulate as-built, pre-damaged and repaired
prestressed concrete girders with the FEM. Kannel et al. 1997 and Okumus et al. 2012 employed the FEM to simulate crack propagation at the ends of girders during prestressing. However, these studies are mainly focused on the transfer of prestressing force to the concrete and ignore the presence of prestressing strands or service loads. O. Yapar et al. (2015) presented a slip-bond failure model for prestressed strands. In this section, we first simulate GDOT girders with O. Yapar et al.'s method. We further simulate pre-damaged and repaired-PSC girders, to understand the effect of reparation on PSC loading capacity.
4.2.1 Damage propagation during pre-stressing and loading We simulate the small-scale beam tests conducted by Yapar et al. (2015). The beam
had a 2590.8 mm length, a 254 mm height and a 99 kN load capacity, as shown in Figure 76. Two strands 12.70 mm in diameter and with 1862 MPa ultimate tensile strength were placed at the bottom of the beam. The top rebar and stirrups had a diameter of 9.53 mm (#3) and 5.08 mm (#2), respectively. The beam was tested under four-point bending condition, as shown in Figure 77. Considering the symmetry of the problem, we model half of the beam length.
117

Figure 76. Details of prestressed concrete test beam (mm) (O. Yapar et al. 2015)
Printed using Abaqus/CAE on: Fri Aug 21 13:43:10 Eastern Daylight Time 2020
Symmetry Figure 77. Applied loading configuration of the test beam (mm)
118

To simulate realistic stress transfer from strands to concrete, an exponential traction-separation behavior is adopted to model the concrete-steel interface, as shown in Figure 78. The bond parameters of the strands and rebars are taken from the study of O. Yapar et al. Parameters are reported in Tables 10 and 11.

Figure 78. Exponential traction separation law (mm) (O. Yapar et al. 2015)

Table 10. Bond parameters of the strands.

Parameter Stiffness of interface Friction coefficient
Normal behavior Maximum separation for damage initiation
Total displacement Damage evolution exponential parameter

Value 47,884 N/mm
0.4 Hard contact
0.42 mm 25.4 mm
4.3

119

Table 11. Bond parameters of the reinforcing rods.

Parameter Stiffness of interface Friction coefficient
Normal behavior Maximum separation for damage initiation
Total displacement Damage evolution exponential parameter

Value 19,064 N/mm
0.4 Hard contact
0.42 mm 25.4 mm
4.3

The beam is modeled with CDP elements, with the material parameters calibrated in Chapter 2. The strands are modeled with 3D solid elements (instead of truss elements in the RC beam model presented in Chapter 3) to achieve better accuracy when simulating the stress release between concrete and the strands. The simulations are performed in three steps, as follows:
Step 1: The bottom two strands are stretched (tensioned) while the bond between the strands and the concrete are deactivated (0.8Y- = 1325) as shown in Figure 79 (a). Step 2: The interfacial bond between the strand and the concrete is restored, and the strands are released to simulate prestress transfer from the strands to the concrete. Step 3: The external force is applied to the beam to the four-point bending test. The details of the concrete and steel meshes are shown in Figure 79. Concrete was modeled with C3D4 linear tetrahedron elements with an edge size of 10 mm. From Figure 79 (a), we can see that, as expected, concrete does not undergo any stress in Step 1. During the prestressing transfer in Step 2, damage develops at the ends of the girder, as shown in Figure 79 (b). Figures 76 (c) and (d) show the stress state induced by the release of the prestressing force in the concrete and in the steel, respectively.

120

(a)
(b)
(c)
(d) Figure 79. Simulation result of prestressed concrete girder (a) stretched strand to initial prestressing force level (b) longitudinal endzone damages followed by releasing of the
strands (c) stress behavior of concrete followed by releasing of the strands (d) stress behavior of rebar and strands followed by releasing of the strands.
121

Figure 80 shows the results of the four-point bending test on pretensioned prestressed concrete girders. The load-displacement curve is in good agreement with that obtained in the experiments reported in (O. Yapar et al. 2015). To assess the effect of prestressing on girder performance, we simulated the four-point bending test without prestressing force. Figure 80(b) shows that the prestressed girder has a loading capacity of 130 kN, while the girder without prestressing force withstand roughly 80 kN. The stress distributions (shear and effective stress) along the length of the strand at the end of girder are shown in Figures 80 (c) and (d). The transferred shear stress from the strands to the concrete exhibits a peak shear stress at a distance of roughly 100 mm from the girder end, which is followed by a plateau after the shear stress reaches a range of -10 to 0 MPa. The curve is roughly hyperbolic. The von Mises stress, defined as & : (with s being the
D
deviatoric stress), reaches 1200 MPa at a distance 380 mm from the girder end. To summarize, Yapar's method can be used to successfully capture the propagation of damage in PSC beams, both during pre-stressing and loading. We use this method in the following to assess the effect of pre-damage and reparation on PSC beam performance.
122

Force [kN]

150

150

Force [kN]

100

100

50

50

FE-PSC

FE-PSC

Experiment

NO-Prestress

0

0

0

5

10

15

20

0

5

10

15

20

Displacement [mm]

Displacement [mm]

(a)

(b)

30

FE-Strand-shear

1200

1000
20 800

Stress [MPa]

10

0

-10 0

200 400 600 800 1000 1200

Distance from girder end [mm]

(c)

600

400

200
0 0

FE-Strand
200 400 600 800 1000 1200
Distance from girder end [mm]
(d)

Figure 80. Simulation of a concrete girder subjected to prestressing followed by fourpoint bending. (a) Comparison of the force vs displacement curves obtained numerically to that reported in Yapar's experiments. (b) Force vs displacement curves with prestress
and without prestress (FEM results). (c) Transferred shear stress from strands to concrete. (d) Effective stress in the strands.

Stress [MPa]

123

4.2.2 Simulation of pre-damaged PSC girders We generate initial damage in the girder model presented in the previous section,
to reproduce four types of end zone damage described in the NCHRP and the PCI guidelines, as shown in Figure 81. In each of the four cases, the total volume of damaged concrete is 75,000 mm3 and the crack width is 5 mm. In the damaged zone, the tensile strength of concrete is assigned 10% of the GDOT tensile strength (i.e., 10% of 2.28 MPa). We subject the pre-damaged girders to pre-stressing, stress release followed by four-point bending. The final damage distributions after Step 3 are shown in Figure 82. Figure 83 compares the load-displacement curves of the four pre-damaged girders in reference to the as-built PSC girder. The four pre-damaged PSC beams exhibit a loss of stiffness of 20 to 25%. The load capacity is almost identical to the load capacity of the as-built PSC beam. The highest loss of capacity occurs when longitudinal cracks exist at the end of the PSC girders.
Figure 81. End zone damage configurations described by the PCI (a) Diagonal crack dia (b) Flange bottom - bot (c) Vertical crack - ver (d) Longitudinal crack - longi
124

(a)

(b)

(c)

(d)

Figure 82. Damage distribution after pre-stressing, strand release and four point

bending of PSC girders with initial damaged end zones. (a) Diagonal crack - dia (b)

Flange bottom crack - bot (c) Vertical crack - ver (d) Longitudinal crack longi

150

Force [kN]

100

50

FE-PSC

Damage-bot

Damage-dig

Damage-longi

Damage-ver

0

0

5

10

15

20

Displacement [mm]

Figure 83. Load-displacement curves of PSC girders without and with initial damage.

125

4.2.3 Simulation of repaired PSC girders with pre-stress damage We use the PSC model presented in section 4.2.1 to calculate the field of maximum
principal strain after applying the pre-stress and determine repairable concrete zone. Figure 84 (a) shows that prestressing causes the creation of longitudinal and diagonal damage zones. Inside these zones, we extract the area where the crack width is at least 0.01 mm, which is the minimum width necessary to allow HMWM penetration, according to our RC beam simulations presented in Section 4.1. The area in which cracks have a width of at least 0.01 mm is represented in grey color in Figure 84 (b). We use MATLAB to binarize the image representing the crack zones (Figure 84 (c)). The white and black colors are identified as damaged and non-damaged areas by the tracing image tool in SOLIDWORKS (Figure 84 (d)). We then simulate the four-point bending test with the PSC beam that has been damaged by prestressing and then repaired. Figure 85 shows that if the cracks induced by pre-stressing are repaired by HMWM deposition, then the PSC girder has a loading capacity that is about 7% higher than that of the as-build PSC girder. Lastly, we simulate the four-point bending with a RC girder that has the same volume fraction of steel rebar than the PSC girder under study, in the case when the RC girder is initially non-damaged, and in the case when the RC girder has a repaired initial damage zone similar to the one that developed during prestress in the PSC girder. Figure 86 shows that the load capacity of the as-built RC girder is about half of that of the as built PSC girder, and that of the damaged-and-repaired RC girder is 75% of that of the as built PSC girder (or 150% of that of the as built RC girder). This comparison shows that HMWM deposition has the potential to strengthen both PSC and RC beams, and that reparation is particularly effective in
126

members with low prestress or low rebar volume fraction (as was the case in the RC beams tested in the laboratory and studied in Chapter 3).

(a) (b) (c)

Crack width 0.01mm (Cracks at 80kN)

(d) Figure 84. Post-processing of simulation result on PSC girder (a) Damage after releasing prestressed strands (b) Zone where crack width is at least 0.01 mm after
releasing prestressed strands (grey color) (c) Binary image using a 0.01mm crack width threshold (d) Concrete part with cracks less than 0.01 mm in width.

127

150

Force [kN]

100

50

FE-PSC

Repaired with Prestressing force

0

0

5

10

15

20

Displacement [mm]

Figure 85. Load vs. displacement curve of as-built PSC girder and PSC girder with initial pre-stressed damaged zone repaired by HMWM deposition

150

Force [kN]

100

50

FE-PSC

Repaired without Prestressing force

Nonrepaired without Prestressing force

0

0

5

10

15

20

Displacement [mm]

Figure 86. Load vs. displacement curve of as-built PSC girder, RC girder of same rebar volume fraction and RC girder with initial pre-stressed damaged zone repaired by HMWM deposition

128

4.2.4 Simulation of repaired PSC girders with end-zone damage We now compare the load capacity of damaged and repaired PSC girders to that of
the as-built PSC girders. Here, we consider the four damage end zones studied in in Section 4.2.2. We replace the damaged elements by repaired CDP elements in all the areas where the crack width exceeds 0.01 mm. The load-displacement curves obtained after simulating pre-stressing, strand release and four-point bending loading are shown in Figure 87. We can see that all pre-damaged-and-repaired beams recovered some stiffness and load capacity. The loss of stiffness in the pre-damaged beams was 20 to 25%. After repair, the loss of stiffness is less than 5%, and the loading capacity is greater than or equal to the asbuilt loading capacity of the as-built girder. We conclude that HMWM deposition allows recovery or mechanical properties in PSC girders that exhibit pre-stress damage or end zone damage as described by the PCI.
150

Force [kN]

100

50

FE-PSC

Repaired-bot

Repaired-dig

Repaired-longi

Repaired-ver

0

0

5

10

15

20

Displacement [mm]

Figure 87. Load-displacement curves of PSC girders without and with initial damage repaired by HMWM deposition.

129

CHAPTER 5. SUMMARY AND CONCLUSIONS
Phase II of this research projects aimed to simulate fracture patterns in concrete, and to numerically assess the effect of HMWM deposition on mechanical properties of concrete structural members.
We calibrated two numerical models to simulate fracture patterns in concrete: a DEM model for analyzing crack patterns at the aggregate scale, and the CDP model for predicting fracture propagation at the laboratory specimen scale. The main conclusions of this part of the project are the following:
1. In both the splitting test and in the uniaxial compression test, a linear relationship exists between the contact ratio a in the ITZ and concrete strength, for 0.5 < a < 1.0.
2. Concrete strength is linearly related to aggregate tensile strength. Weak aggregates decrease concrete strength. The aggregate volume fraction is usually around 30% or more. We found that multiplying the aggregate tensile strength by four could increase concrete strength by 2/3. However, the influence of aggregate tensile strength on concrete strength is much lower than that of the ITZ.
3. The CDP model can predict fracture patterns that develop in concrete in tension and compression. This model is used for simulating repaired concrete tests.
We then studied concrete/epoxy strength at molecular and metric scales. The main conclusions are the following:
4. We built MD concrete/epoxy and concrete/HMWM interface models by constructing a layer of crystalline silica and calcite adjacent to a layer of highly cross-linked epoxy. We simulated cross-linking of DGEBA with DETA and cross-linking of HMWM. We then simulated annealing and NPT steady state, which allowed us to build an MD epoxy model
130

with the desired thermal expansion coefficient, transition temperature to glassy state and mass density. MD simulations of tension tests allowed us to determine the elastic properties of the epoxy models. It was found that cross-linked epoxy mechanical properties highly depend on the conversion ratio (i.e. on the cross-linking rate), the specimen size and the loading rate.
5. MD simulations of uniaxial tension tests showed that all interfaces fail by debonding, at the surface of contact between HMWM and the mineral substrate, and that the interfacial strength decreases in the presence of moisture, under low strain rate, or at high temperature. Silica/HMWM interfaces were systematically stronger than the calcite/HMWM interfaces. The results showed that the work of separation was mostly done by van der Waals forces.
6. We showed our MD results could be extrapolated to predict the tensile modulus of HMWM at low strain rate. The Johnson-Cook model was used to predict the HMWM/mineral interfacial strength for a broad range of temperatures and strain rates. MD results confirm that HMWM should be applied on dry surfaces and in concrete exposed to mild temperatures. HMWM application in concrete with high silica contents is more likely to last than in concrete with high calcite contents.
7. We calibrated material properties of repaired concrete with HMWM against repaired uniaxial compression tests and three-point bending tests. We assigned the calibrated CDP model to concrete volume elements and concrete/HMWM interfaces were modelled with CZ elements with a bilinear softening law. The FEM model enriched with CZ elements was used to simulate repaired reinforced concrete members.
The last part of the study focused on RC and PSC beams. The main findings are as follows:
131

8. We used DIC to identify the zones of high maximum principal strain after the first loading cycle applied to RC beams tested in Phase III. These zones were assigned repaired concrete elements and HMWM CZ elements for simulating the second load cycle, after reparation. The distribution of damage and the field of maximum principal strain calculated numerically at failure match the damage zone and the maximum principal strain field found experimentally. Simulation results suggest that HMWM can penetrate cracks of width 0.01 mm and above by gravity. We also found that HMWM reparation increase concrete stiffness and strength if crack in concrete members are over 0.1 mm in width, in which case, the load capacity of repaired RC beams is 30 to 40% higher than that of as-built RC beams.
9. We simulated pre-stressing, strand release and four-point loading of PSC girders. The FEM model matched results reported in the literature and captured the development of damage during pre-stressing. We found that the load capacity of a PSC girder damaged by pre-stressing and then repaired would be about 7% higher than that of the as-built PSC girder. At same volume fraction of rebar, a RC girder has a load capacity that is about half of that of the PSC girder. Addition of HMWM in a zone of damage similar to the prestressing damage zone would increase the load capacity of the RC girder by 50%, reaching 75% of the load capacity of the as-built PSC girder. We simulated pre-stressing followed by strand release and four-point bending in PSC girders with initial end zone damaged described in the PCI guidelines. Compared to the as-built girders, pre-damaged girders exhibited a loss of stiffness of 20 to 25 % and a loss of load capacity of up to 15%.
132

Reparation allowed recovery of both stiffness and strength. Repaired girders sometimes exhibited a higher load capacity than the as-built PSC girders.
Based on the conclusions drawn in this study, we provide the following design recommendations, supported by the strength charts in Fig. 88-91:
1. Improve concrete Use higher aggregate strength Use smaller aggregates Increase silica contents
2. Optimize HMWM deposition technique Apply at low temperature Apply on dry crack faces Repair cracks > 0.01 mm in width
3. Do repair PSC beams according to PCI guidance Mechanical properties of PSC beams can be recovered upon reparation of cracks classified by the PCI Reparation has a greater influence on load-displacement curve for lower rebar volume fraction and lower prestress Load capacity of repaired RC beams with no prestress can exceed that of as-built RC beams
4. Plot strength charts for in-situ implementation Full scale testing is needed to validate repaired PSC beam simulations
133

Full scale parametric studies is needed to assess a larger set of crack distributions and reinforcement/pre-stress designs
Figure 88. Effect of aggregate strength (top row) and interface transition zone contact ratio (bottom row) on concrete tensile strength (left) and compressive strength (right)
134

Load [kN] Load [kN]

Figure 89. Effect of concrete mineral composition, temperature and moisture on concrete/HMWM interface strength.

300

Repaired after major crack #1

250

Repaired after major crack #2 Crack repaired > 0.01mm

Crack repaired > 0.05mm

200

Crack repaired > 0.1mm

300

RC concrete #1

250

RC concrete #2 Repaired after minor crack

Crack repaired > 0.01mm

200

Crack repaired > 0.005mm

150

150

100

100

50

50

0

0

0

5

10

15

0

5

10

15

Displacement [mm]

Displacement [mm]

Figure 90. Effect of HMWM reparation on the load capacity of RC beams

135

Force [kN] Force [kN]

150

150

100

100

50

FE-PSC

50

FE-PSC

Damage-bot

Repaired-bot

Damage-dig

Repaired-dig

Damage-longi

Repaired-longi

Damage-ver

Repaired-ver

0

0

0

5

10

15

20

0

5

10

15

20

Displacement [mm]

Displacement [mm]

Figure 91. Effect of HMWM reparation on the load capacity of PSC beams

136

CHAPTER 6. REFERENCE
ACIMP - American Concrete Institute Manual of Practices. Crack widths in concrete bridge deck. ACI 1997;224R:4.4.
Ahmad, S., & Alghamdi, S. A. (2012). A Study on Effect of Coarse Aggregate Type on Concrete Performance. Arabian Journal for Science and Engineering, 37(7), 1777-1786. doi:10.1007/s13369-012-0282-6 Aitcin, P. C., & Mehta, P. K. (1990). Effect of Coarse-Aggregate Characteristics on
Mechanical-Properties of High-Strength Concrete. ACI Materials Journal, 87(2), 103-107. Alder, B. J., & Wainwright, T. E. (1959). Studies in molecular dynamics. General method. Journal of chemical physics, 31(2), 459-459. doi:10.1063/1.1730376 Arab, B., Shokuhfar, A., & Ebrahimi-Nejad, S. (2012). V. Proceedings of the International Conference Nanomaterials: Applications and Properties, 1(1), 1-4. ASTM, C. C39M (2012) Standard test method for compressive strength of cylindrical concrete specimens. American Society for Testing and Materials, ASTM, West Conshohocken.
ASTM. A706/A706M - 15 Standard Specification for Low-Alloy Steel Deformed and Plain Bars for Concrete. Astm Int 2015:16. https://doi.org/10.1520/A0706_A0706M-16.
ASTM. Standard Specification for Deformed and Plain Carbon-Steel Bars for Concrete Reinforcement. A615/A615M - 18 2016:18. https://doi.org/10.1520/A0615.
Attanayake U, Xuemei L, Simon Ng PE, Haluk Aktan, P.E. MA. Penetrating Sealants for Concrete Bridge Decks--Selection Procedure. J Bridg Eng 2006;11(5):53340. https://doi.org/10.1061/(ASCE)1084-0702(2006)11.
Bandyopadhyay, A., Valavala, P. K., & Odegard, G. M. (2009). Atomistic Modelling of Crosslinked Epoxy Polymer. Structure(17th ICCM Conference Proceedings).
Bardella, L. (2001). A phenomenological constitutive law for the nonlinear viscoelastic behaviour of epoxy resins in the glassy state. European Journal of Mechanics aSolids, 20(6), 907-924. doi:Doi 10.1016/S0997-7538(01)01180-9
Basunbul IA, Gubati AA, Al-Sulaimani GJ, Baluch MH. Repaired reinforced concrete beams. ACI Mater J 1990;87:34854. https://doi.org/10.14359/2064.
Bazant ZP, Becq-Giraudon E. Statistical prediction of fracture parameters of concrete and implications for choice of testing standard. Cem Concr Res 2002;32:52956. https://doi.org/10.1016/S0008-8846(01)00723-2.
Bazant, Z. P., & Pfeiffer, P. A. (1987). Determination of Fracture Energy from Size Effect and Brittleness Number. ACI Materials Journal, 84(6), 463-480.
137

Bazant, Z. P., Kazemi, M. T., Hasegawa, T., & Mazars, J. (1991). Size effect in Brazilian split-cylinder tests: measurements and fracture analysis. ACI Materials Journal, 88(3), 325-332.
Beushausen, H., & Dittmer, T. (2015). The influence of aggregate type on the strength and elastic modulus of high strength concrete. Construction and Building Materials, 74, 132-139. doi:10.1016/j.conbuildmat.2014.08.055
Biel TD, Lee H. Performance study of Portland cement concrete pavement joint sealants. J Transp Eng 1997;123:398404. https://doi.org/10.1061/(ASCE)0733947X(1997)123:5(398).
Blaber J, Adair B, Antoniou A. Ncorr: Open-Source 2D Digital Image Correlation Matlab Software. Exp Mech 2015;55:110522. https://doi.org/10.1007/s11340015-0009-1.
Camanho P, Davila CG. Mixed-Mode Decohesion Finite Elements in for the Simulation Composite of Delamination Materials. 2002. https://doi.org/10.1177/002199803034505.
Chi, J. M., Huang, R., Yang, C. C., & Chang, J. J. (2003). Effect of aggregate properties on the strength and stiffness of lightweight concrete. Cement & Concrete Composites, 25(2), 197-205. doi:Pii S0958-9465(02)00020-3
Committee, A. (1999). Building code requirements for structural concrete:(ACI 318-99); and commentary (ACI 318R-99).
Coronado CA, Lopez MA. Experimental characterization of concrete-epoxy interfaces. J Mater Civ Eng 2008;20:30312. https://doi.org/10.1061/(ASCE)08991561(2008)20:4(303).
Cundall, P. A. (1971). A computer model for simulating progressive, large-scale movements in blocky rock systems. Paper presented at the Proc. Int. Symp. on Rock Fracture.
Cundall, P. A., & Strack, O. D. L. (1979). Discrete Numerical-Model for Granular Assemblies. Geotechnique, 29(1), 47-65.
Ding, X. B., & Zhang, L. Y. (2014). A new contact model to improve the simulated ratio of unconfined compressive strength to tensile strength in bonded particle models. International Journal of Rock Mechanics and Mining Sciences, 69, 111-119. doi:10.1016/j.ijrmms.2014.03.008
Doi 10.1016/S0958-9465(02)00020-3
Dry C, Corsaw M, Bayer E. A comparison of internal self-repair with resin injection in repair of concrete. J Adhes Sci Technol 2003;17:7989. https://doi.org/10.1163/15685610360472457.
138

Dry CM. Repair and prevention of damage due to transverse shrinkage cracks in bridge decks. Smart Struct Mater 1999 Smart Syst Bridg Struct Highw 1999;3671:2536. https://doi.org/10.1117/12.348675.
Dry CM. Three designs for the internal release of sealants, adhesives, and waterproofing chemicals into concrete to reduce permeability. Cem Concr Res 2000;30:196977. https://doi.org/10.1016/S0008-8846(00)00415-4.
Edward M. Wu R. C. Router J. Crack Extension in Fiberglass Reinforced Plastics. Univ Illinois 1965;T. & AM Re.
Elices M, Guinea G V., Gmez J, Planas J. The cohesive zone model: Advantages, limitations and challenges. Eng Fract Mech 2001;69:13763. https://doi.org/10.1016/S0013-7944(01)00083-2.
Fakhimi, A., & Villegas, T. (2007). Application of dimensional analysis in calibration of a discrete element model for rock deformation and fracture. Rock Mechanics and Rock Engineering, 40(2), 193-211. doi:10.1007/s00603-006-0095-6
Gallucci, E., Scrivener, K., Groso, A., Stampanoni, M., & Margaritondo, G. (2007). 3D experimental investigation of the microstructure of cement pastes using synchrotron X-ray microtomography (mu CT). Cement and Concrete Research, 37(3), 360-368. doi:10.1016/j.cemconres.2006.10.012
Glvez, J. C., Cervenka, J., Cendn, D. A., & Saouma, V. (2002). A discrete crack approach to normal/shear cracking of concrete. Cement and Concrete Research, 32(10), 1567-1585. doi:10.1016/S0008-8846(02)00825-6
Giaccio, G., Rocco, C., & Zerbino, R. (1993). The Fracture Energy (G(F)) of High-Strength Concretes. Materials and Structures, 26(161), 381-386. doi:Doi 10.1007/Bf02472938
Hayes AR. The Air Force Institute of Technology. IRE Trans Educ 1962;5:1178. https://doi.org/10.1109/TE.1962.4322266.
Hillerborg, A., Moder, M., & Petersson, P.-E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 6(6), 773-781.
Issa, C. A., & Debs, P. (2007). Experimental study of epoxy repairing of cracks in concrete. Construction and Building Materials, 21(1), 157-163. doi:10.1016/j.conbuildmat.2005.06.030
Itasca. (2008). PFC3D manual, version 4.0. In: Itasca Minneapolis.
Jankowiak T, Lodygowski T. Identification of parameters of concrete damage plasticity constitutive model. Found Civ Environ ... 2005:5369.
139

Jin W, Xu H, Arson C, Busetti S. Computational model coupling mode II discrete fracture propagation with continuum damage zone evolution. Int J Numer Anal Methods Geomech 2017;41:22350. https://doi.org/10.1002/nag.2553.

Labibzadeh M, Zakeri M, Shoaib AA. A new method for CDP input parameter identification of the ABAQUS software guaranteeing uniqueness and precision. Int J Struct Integr 2017;8:26484. https://doi.org/10.1108/IJSI-03-2016-0010.

Landis, E. N., & Nagy, E. N. (2000). Three-dimensional work of fracture for mortar in compression. Engineering Fracture Mechanics, 65(2-3), 223-234. doi:Doi 10.1016/S0013-7944(99)00124-1
Lau, D., Bykztrk, O., & Buehler, M. J. (2012). Characterization of the intrinsic strength between epoxy and silica using a multiscale approach. Journal of Materials Research, 27(14), 1787-1796. doi:10.1557/jmr.2012.96

Lee J, Fenves GL. Plastic-damage model for cyclic loading of concrete structures. J Eng Mech 1998;124:892900. https://doi.org/10.1061/(ASCE)07339399(1998)124:8(892).

Li, C., & Strachan, A. (2010). Molecular simulations of crosslinking process of

thermosetting

polymers.

Polymer,

51(25),

6058-6070.

doi:10.1016/j.polymer.2010.10.033

Li, D. Y., & Wong, L. N. Y. (2013). The Brazilian Disc Test for Rock Mechanics Applications: Review and New Insights. Rock Mechanics and Rock Engineering, 46(2), 269-287. doi:Doi 10.1007/S00603-012-0257-7

LIANG Y-C. Experimental study and theoretical modeling of recycled aggregate concrete and evaluation of long-term performance of reinforced concrete bridge decks. 2012.

Lubliner J, Oliver J, Oller S, Oate E. A plastic-damage model for concrete. Int J Solids Struct 1989;25:299326. https://doi.org/10.1016/0020-7683(89)90050-4.

Ma, Y. F., & Huang, H. Y. (2018). A displacement-softening contact model for discrete element modeling of quasi-brittle materials. International Journal of Rock Mechanics and Mining Sciences, 104, 9-19. doi:10.1016/j.ijrmms.2018.02.007

Mander JB, Priestley MJ, Park R. Theoretical stress-strain model for confined concrete. J Struct Eng (United States) 1988;114:180426. https://doi.org/10.1061/(ASCE)0733-9445(1988)114:8(1804).

Martin J, Stanton J, Mitra N, Lowes LN. Experimental testing to determine concrete fracture energy using simple laboratory test setup. ACI Mater J 2007;104:57584. https://doi.org/10.14359/18961.

140

Meggers DA. Final Report Crack Sealing and Repair of Older Serviceable Bridges Using Polymer Sealers. 1998. https://doi.org/FHWA-KS-98-4.
Michael M. REPAIR AND PROTECTION OF HYDRAULIC CEMENT CONCRETE BRIDGE DECKS Sprinkel Manager Cristopher Standard Title Transportation. 1994.
Park K, Ha K, Choi H, Lee C. Prediction of interfacial fracture between concrete and fiber reinforced polymer (FRP) by using cohesive zone modeling. Cem Concr Compos 2015;63:12231. https://doi.org/10.1016/j.cemconcomp.2015.07.008.
Parks, G. A. (1984). Surface and Interfacial Free-Energies of Quartz. Journal of Geophysical Research, 89(Nb6), 3997-4008. doi:DOI 10.1029/JB089iB06p03997
PETERSSON P-E. CRACK GROWTH AND DEVELOP MENT OF FRACTURE ZONES IN PLAIN CONCRETE AND SIMILAR MATERIALS. Sweden: 1981. https://doi.org/10.5144/0256-4947.1998.84.
Petersson, P. E. (1980). Fracture Energy of Concrete - Method of Determination. Cement and Concrete Research, 10(1), 78-89.
Potyondy, D. O., & Cundall, P. A. (2004). A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, 41(8), 1329-1364.
Rahim A, Jansen D, Abo-Shadi N, Simek J. Overview of high-molecular-weight methacrylate for sealing cracks in concrete bridge decks. Transp Res Rec 2010;2202:7781. https://doi.org/10.3141/2202-10.
Reneson KJAESCEFJ. Crack and Concrete Deck Sealant Performance Crack and Concrete Deck Sealant Performance. 2009.
Rocco C, Guinea G V., Planas J, Elices M. Review of the splitting-test standards from a fracture mechanics point of view. Cem Concr Res 2001;31:7382. https://doi.org/10.1016/S0008-8846(00)00425-7.
Rocco, C., Guinea, G. V., Planas, J., & Elices, M. (1999). Size effect and boundary conditions in the Brazilian test: Experimental verification. Materials and Structures, 32(217), 210-217. doi:Doi 10.1007/Bf02481517
Sasse, H. R., & Fiebrich, M. (1983). Bonding of polymer materials to concrete. Mat??riaux et Constructions, 16(4), 293-301. doi:10.1007/BF02473695
Schreppers GJ, Frissen C, Kang HJ. Prediction of crack-width and crack-pattern. TNO Diana 2011.
Scrivener, K. L., Crumbie, A. K., & Laugesen, P. (2004). The interfacial transition zone (ITZ) between cement paste and aggregate in concrete. Interface Science, 12(4), 411-421. doi:DOI 10.1023/B:INTS.0000042339.92990.4c
141

Seok S, Haikal G, Ramirez JA, Lowes LN. High-resolution finite element modeling for bond in high-strength concrete beam. Eng Struct 2018;173:91832. https://doi.org/10.1016/j.engstruct.2018.06.068.
Simulia DS. Abaqus 6.12 documentation. Rhode Island: 2012.
Skarzynski, L., & Tejchman, J. (2016). Experimental Investigations of Fracture Process in Concrete by Means of X-ray Micro-computed Tomography. Strain, 52(1), 26-45. doi:10.1111/str.12168
Snoeck D, Van Tittelboom K, De Belie N, Steuperaert S, Dubruel P. The use of superabsorbent polymers as a crack sealing and crack healing mechanism in cementitious materials. Concr Repair, Rehabil Retrofit III - Proc 3rd Int Conf Concr Repair, Rehabil Retrofit ICCRRR 2012 2012:1527.
Sprinkel MM, Sellars AR, Weyers RE. Rapid Concrete Bridge Deck Protection, Repair and Rehabilitation, SHRP-S-344, Concrete Bridge Protection and Rehabilitation: Chemical and Physical Techniques. 1993.
State of Georgia D of T. Bridge and Structures Design Manual. Bridg Struct Des Man 2015.
Transpo Industries. Sealate T-70 Technical Data Sheet. 2017.
Truster TJ. Discontinuous Element Insertion Algorithm. Fac Publ Other Work -- Civ Environ Eng 2015.
Wayne D. Mangum, Alejandro J. Bermudez-Goldman, David P. Whitney, David W. Fowler and AHM. REPAIRING CRACKS IN PORTLAND CEMENT CONCRETE USING POLYMERS. vol. 48. 1986. https://doi.org/10.1016/0041008X(79)90032-2.
Weyers R, Prowell B, Sprinkel M, Vorster M. Concrete bridge protection, repair, and rehabilitation relative to reinforcement corrosion: A methods application manual. vol. 100. Strategic. National Academy of Sciences; 1993.
Wohl RL, LaFraugh RW. Criteria for the selection of penetrating hydrophobic sealers used in the repair of concrete parking decks. ASTM Spec Tech Publ 1990:7582. https://doi.org/10.1520/stp26089s.
Xiong, G., Liu, J., Li, G., & Xie, H. (2002). A way for improving interfacial transition zone between concrete substrate and repair materials. Cement and Concrete Research, 32(12), 1877-1881. doi:10.1016/S0008-8846(02)00840-2
Xu, Y., Xu, C. H., Zhou, Z., Du, J., & Hu, D. P. (2010). 2D DEM simulation of particle mixing in rotating drum: A parametric study. Particuology, 8(2), 141-149. doi:10.1016/j.partic.2009.10.003
142

Yang, S., Gao, F., & Qu, J. (2013). A molecular dynamics study of tensile strength between a highly-crosslinked epoxy molding compound and a copper substrate. Polymer (United Kingdom), 54(18), 5064-5074. doi:10.1016/j.polymer.2013.07.019
Yapar O, Basu PK, Nordendale N. Accurate finite element modeling of pretensioned prestressed concrete beams. Eng Struct 2015;101:16378. https://doi.org/10.1016/j.engstruct.2015.07.018.
Zhou, F. P., Lydon, F. D., & Barr, B. I. G. (1995). Effect of Coarse Aggregate on ElasticModulus and Compressive Strength of High-Performance Concrete. Cement and Concrete Research, 25(1), 177-186. doi:Doi 10.1016/0008-8846(94)00125-I
143