Mechanical Properties and Microstructure of Reinforced Passivation Film: A Reactive Molecular Dynamics Study

Reinforced concrete is one of the most critical composite materials in the modern civil engineering and can improve the tensile resistance of concrete. Its passivation film plays an important role in the durability of concrete and the steel corrosion. But, due to the size limitations, the destruction of micro-scale steel bars has not been well studied. In this work, the reactive molecular dynamics simulation was employed to studying the mechanical properties of the steel and its passivation film. The uniaxial stretching of different compounds of γ-FeOOH, γ-Fe2O3 and Fe was performed. We found that the oxidation can reduce the tensile strength of steel. For the three compounds of γ-FeOOH, γ-Fe2O3 and Fe, the order of tensile strength from high to low is Fe > γ-Fe2O3 > γ-FeOOH. But, the ductility of γ-FeOOH under x direction is increased. The detail microstructure analysis shown that the difference of tensile strength is origin from the coordination in the materials. The two kinds of stretching processes of whole system stretching (in the Fe phase and x direction of γ-FeOOH phase) and partly area stretching (in the Fe2O3 phase and z direction of γ-FeOOH phase) were clarified. The external force is dispersed in whole system stretching but opposite in partly area stretching. This investigation leads to possible new direction for studying the tensile strength of materials, and the strategy of evaluating materials tensile strength can supply valuable information in evaluating and improving the mechanical properties of reinforced concrete.


INTRODUCTION
Reinforced concrete can effectively improve the tensile strength of ordinary Portland concrete (Alexander, 1959;Song and Hwang, 2004;Song et al., 2005;Park et al., 2012), which is widely used in high-rise buildings, bridges and other construction projects (Horvath, 2004;Scrivener and Kirkpatrick, 2008;Anandavalli, 2012;Allen and Iano, 2019). But the rebar failure caused by cyclic loading (Psycharis et al., 2018) great limits the application of reinforced concrete. The studies of mechanical properties of rebar have attracted extensive attention (Kanaev et al., 2018;Bandyopadhyay et al., 2019). In general, the rebar contains two components which are the steel phase and its passivation film. In the alkaline environment of concrete, the passivation film is formed and locates on the surface of the steel (Luo et al., 2015). Reinforced passivation film plays an important role for improving the durability of concrete and preventing corrosion. Therefore, the study on the mechanical properties of steel and its' passivation film is facilitated to further using and developing the reinforced concrete.
Many researches effort to study the mechanical properties of steel bars at the macroscopic scale (Tuutti, 1982;MacGregor et al., 1997;Angst et al., 2017;Gong and Maekawa, 2018;Gong and Maekawa, 2019). Thomas and Ramaswamy (2007) have studied the mechanical properties of steel fiber-reinforced concrete and proposed a model that predicted the experimental test data quite accurately. Further, the mechanical properties for different effects of temperature (Dügenci et al., 2015), craft (Seo et al., 2016) and strength (Iqbal et al., 2015) fiber-reinforced concrete were investigated. These studies have played a crucial role for understanding the relationship among the mechanical properties of steel bars, structure, and environment. However, due to the size limitations, the destructions of steel bars in microscale have not been well clarified (Zhang et al., 2016;Boscato et al., 2018;Saha et al., 2018;Nykyforchyn et al., 2019;Sienko et al., 2019). The fiber in reinforced concrete can be regarded as a layered structure in which the surface is composed by γ-FeOOH, γ-Fe 2 O 3, and iron from outside to inside. The iron phase is covered by the passivation film. The mechanical properties of the passivation film can significantly affect the mechanical properties of fiber. The investigation of mechanical properties of passivation film needs micro-scale perspective, due to the fine structure of passivation film on steel surface.
Recently, the molecular dynamics (MD) simulation is considered as an effective means to study the molecular behavior in nanoscale (Abdolhosseini Qomi et al., 2012;Hou et al., 2018;Wang et al., 2019;Fogg et al., 2019;Wang et al., 2020a;Wang et al., 2020b;Wang et al., 2020c), which can also reveal the micro-structural changes of materials under force loading (Hou et al., 2014;Hou et al., 2015;Hou, 2019;. There are many precedents for the study of crystal mechanical behavior using MD simulation. Li et al. (2019) study the tension behaviors of nanopolycrystalline Cu-Ta alloys through MD simulations, and the influences of several important factors on the mechanical properties of the materials are considered. Wu et al. (2015) used MD methods to study the effects of vacancies and dislocations on the mechanical properties for interfacial deformation processes of methane hydrates. The results show dislocation-free brittle failure in monocrystalline hydrates and an unexpected crossover from strengthening to weakening in polycrystals. Lin et al. (2017) used parallel MD to simulate the mechanical properties of graphene reinforced polymer nanocomposites with temperature dependent. These researches confirm that MD is an effective way to study the mechanical properties of materials. But, for many metallic materials, the breaking process often involves the breaking of the main valence bond. Many traditional force-fields cannot accurately descript the internal interaction during the breaking process of materials. In the process of steel bar destruction, the molecular bonding changes are very complicated, and the traditional MD force field cannot describe the process correctly. Van Duin et al. (2001) have proposed a useful reactive force field that can accurately descript the atomic interaction (Manzano et al., 2012a;Manzano et al., 2012b;Pitman and Van Duin, 2012;Hou et al., 2017;Soria et al., 2018;Ponomarev et al., 2019;, which can be used in studying the destruction process of steel and its' passivation film. Thus, in this work, the reactive MD method was employed for studying the mechanical properties of steel and its' passivation film. The uniaxial stretching processes of different compounds (γ-FeOOH, γ-Fe 2 O 3 , and iron) were performed. We found that oxidation can reduce the strength of iron, but increased ductility under certain condition. The strength evolutions of micro-scale structure in the destruction process of each compounds were clarified. The detail analysis reveals the relationship between internal bonding and mechanical properties of steel. This study presented can supply valuable information in evaluating and improving the mechanical properties of reinforced concrete.

Models
The structure used for the simulation was obtained by amplifying the initial lattice cell (The corresponding Fe, γ-Fe 2 O 3 , and γ-FeOOH.cif is shown in Supplementary Material or the link: https://doi.org/10.5281/zenodo.3957278. https://doi.org/10.5281/ zenodo.3961152) for the Fe, γ-Fe 2 O 3 , and γ-FeOOH as shown in Figure 1. The box size was controlled to approximately 20 Å in x/ y direction and 50 Å in z direction. The box parameters of the expanded Fe were a × b × c 20.04 Å × 20.04 Å × 48.72 Å, α β c 90°, which is containing 1,666 iron atoms. The box parameters of the expanded γ-Fe 2 O 3 were a × b × c 17.44 Å × 20.14 Å × 54. 88 Å, α β c 90°, which is containing 944 Fe atoms and 1,184 O atoms. The box parameters of the expanded γ-FeOOH were a × b × c 49.28 Å × 25 Å × 19.35 Å，α β c 90°, which is containing 640 Fe atoms and 1,280 O atoms and 640 H atoms. The three types of models employed the same simulation conditions to ensure the accuracy of simulation results. The temperature was maintained at 300 K and used the Nose-Hoover thermostat (Nosé, 1984). The atomic trajectories were calculated by using the Verlet algorithm and the simulation time step was set to 0.25 fs.

Simulation Details
The LAMMPS software (Plimpton, 1995) was employed to perform the MD simulations. The reactive force field force field was employed for Fe, γ-Fe 2 O 3 , and γ-FeOOH model, which is widely used to study the strength process of materials. The force field parameters were modeled by using the work of Aryanpour et al. (2010). And the relevant force field parameters can be free acquire in their supporting information. The force field parameters are able to reproduce chemical reactions and examined by quantum calculations and experimental measurements in Aryanpour's work. Periodic boundary conditions were employed. Firstly, energy minimization was achieved to the transport simulation to eliminate high energy regions during model construction. Secondly, all atoms in the model relaxed 10 ps under the NPT ensemble. Finally, 800 ps stretched process was performed at a constant rate of 0.0008 Å/ps in the z direction, and the stretching rate is not set in the x and y directions. It should be noted that a large number of atomic trajectories will be used for later statistical analysis. The configuration is recorded every 100 steps to ensure the statistical stability of the data analysis.

Stress-Strain Curve
The stress-strain curve can be used to analyze the mechanical behavior of the layered structure during stretching and helps to understand the constitutive relationship between stress and strain. In the process of uniaxial stretching simulation, the mechanical properties of Fe, γ-Fe 2 O 3 and γ-FeOOH were evaluated by calculating the stress-strain curve at the tensile rate of 0.0008 Å/ps and running at 800 ps. As shown in Figure 2, the stress-strain curves of the three types of structures were plotted. Since the γ-FeOOH structure has significant anisotropy, thus it is stretched in both three directions of x, y and z, and the stress-strain relationship in three directions is compared. For Fe and γ-Fe 2 O 3 , as shown in Figure 2A, during the tensile process, stress first increases linearly in the elastic stage and subsequently slowly increases to a maximum of 11 Gpa at a strain of 0.1 Å/Å. As the strain increases, the stress drops to 8 Gpa gradually and has a longer yield behavior, which means that Fe has good toughness in the z direction. Although the γ-Fe 2 O 3 structure cannot be completely destroyed by the stretching in the simulation time, it showed a slightly complete destruction process. For γ-Fe 2 O 3 , at strains up to 0.04 Å/Å, the stress reaches the elastic limit (about 6 Gpa). After the strain reaches 0.04 Å/Å, the material stress drops sharply without obvious yield behavior. Even at a strain of 0.16 Å/Å, the stress does not drop to zero, implying good plasticity along the z direction. As shown in Figure 2B, γ-FeOOH exhibits different stress-strain relationships of tensile loading in the x, y, and z direction. Since the y direction mainly depends on the interlayer non-bonding interaction between the interlayer iron and hydroxyl groups, the intermolecular force is weaker than the covalent bond, so the elastic in the y direction is significant lower than that in the x or z direction. The y-direction layers are separated from each other after that the tensile load is applied along z direction. Such performance of stress-strain curve indicates that the oxidation can reduce the tensile strength of steel. For the three compounds of γ-FeOOH, γ-Fe 2 O 3 , and Fe, the order of tensile strength from high to low is Fe > γ-Fe 2 O 3 > γ-FeOOH. But, the ductility under x direction is increased, which indicates that the surficial γ-FeOOH facilitates to resist bends.
It is observed through the stress-strain curves of γ-FeOOH in three directions that the stress-strain curves in the x and y direction have a complete failure process. When stretched in the x direction and the strain is 0.02 Å/Å, the structure reaches the elastic limit and then enters the yielding stage. During the yield phase, the structure loses its ability to resist deformation. As the simulation progresses, the strain continues to increase. When the And the supercell structures of (D) Fe, (E) γ-Fe 2 O 3 , and (F) γ-FeOOH. The purple, red and white balls represent iron, oxygen, and hydrogen atoms, respectively.
Frontiers in Materials | www.frontiersin.org September 2020 | Volume 7 | Article 523126 strain reaches 0.04 Å/Å, the structure enters the strengthening phase. During this period, because of the dislocations, slips, twins and other line defects caused by stress load, the material has the significantly effects on dislocation strengthening, so the stress increases slightly. When the strain reaches 0.09 Å/Å, the material enters the local necking behavior and defects increase, resulting in a sharp drop in stress. The structure reaches the elastic limit when the tensile strain in the y direction reaches 0.05. Then there is a small stage of the yield behavior and the strengthening behavior. Unlike stretching in the x direction, when stretching in the y direction, there is no significant local necking behavior. This is due to the weak hydrogen bonding effect, which cause the interlayer separation process is weak by the contraction of the bond. For further analyzing the influence of local variation on mechanical properties when the material is subjected to load, we calculated relevant parameters of the radial distribution function (RDF), time correlation function and coordination number.

Microstructural Parameter
The RDF is usually used to evaluate the degree of disorder of the structure, reflecting the atomic aggregation state of the material. The atomic arrangement of initial models are explained by calculating the RDF of Fe, γ-Fe 2 O 3 , and γ-FeOOH before stretching. Fe structure contains only iron atom, and adjacent iron atoms are bonding to surrounding atoms by forming the metal bonds. As shown in Figure 3, there is a multi-peak distribution of the RDF curve of Fe-Fe. The valley is located at 2.44 Å, which determines the number of coordination atoms of the central Fe atom in this range, while the valley at 3.2 Å contains the remaining two atoms around the central Fe atoms that are not bonded to the central Fe atoms. The primary bonding in γ-Fe 2 O 3 structure is the Fe-O s (O s is a bridging oxygen atom in the structure) bond and Fe-Fe bond. In the steady state before stretching, there is a small valley at 1.97 Å and a large valley at 2.6 Å between the RDF curve of Fe-O s . The locations of two valleys are the distances of the bonding length of Fe-O s bond that does not form a regular tetrahedron centered on Fe atoms exhibited the anisotropy. In the structure of γ-FeOOH, Fe-O s and Fe-O h rely on the covalent bond with strong interaction. γ-FeOOH exhibits a distinct layered distribution in the y direction, and the interaction of layers is mainly driven by the intermolecular force: Van der Waals force and H-bond that results the relatively weak interactions.
As shown in Figure 4, the change in atomic and local environment during stretching is revealed by calculating the RDF curve changes corresponding to different stretching stages. In the Fe structure, before the strain is 0.04 Å/Å, the RDF curve of Fe-Fe showed the same tendency as the steady state before stretching: there are two valleys at the distance of 2.44 and 3.2 Å. However, as the strain increases, the depth of the first valley in the curve becomes shallower. When the strain reaches 0.12 Å/Å, the first valley disappears completely. At the same time,  Frontiers in Materials | www.frontiersin.org September 2020 | Volume 7 | Article 523126 the first peak width in the RDF curve of Fe-Fe is gradually widened during the stretching process, indicating that the structure tends to be disordered by the relatively regular crystal during the stretching process. The bond between Fe-Fe at 2.44 Å during stretching is first elongated, while the covalent bond with a bond length of 3.2 Å is more stable. Similarly, in the γ-Fe 2 O 3 structure, the location of first valley of RDF gradually decreases during the process of strecthing. The bond with the bond length of 1.98 Å are gradually elongated, while local environments with 2.5 Å are more stable. As the simulation progresses, the RDF curve is tended to flat and the numbers of peaks is reduced, which means that the structure changes from crystal to amorphous form in the stretching progresses. On the contrary, as the stretching progresses, the RDF curve of Fe-O s in the γ-FeOOH structure does not change much. However, as the width of the first peak gradually increase as the stretching proceeds, meaning that the structure tends to be disordered. In γ-FeOOH, the first peak (1.98 Å) indicates the most stable bond distance between Fe-O s , and the second peak (3.2 Å) indicates the strong spatial correlation between neighboring Fe atoms.
The different RDFs in the stretching process proof that the orderly of all the three materials are decreased after stretching. But the ability of maintaining the initial ordering of Fe phase is better than that of γ-Fe 2 O 3 and γ-FeOOH phase, which indicates the highest tensile strength appeared on Fe phase. The RDF of γ-Fe 2 O 3 phase has the largest change, exhibiting the lowest ductility. We can conclude that the stretching of steel primary affects the structural ordering. But, the first peak of RDF seems little changes, which indicates that the lattice is less affected in the directions where a stretching rate is not set.

Strain and Coordination
For exploring the origin of different tensile strength in three model, Figure 5 shows the change in atomic coordination number of Fe during simulation. The time evolution of microstructure is consistent with the RDF analysis. By analyzing the changes of coordination number, reveal the change of the local environment of the atom and the structural during the stretching process. During the simulation, the atoms inside the crystal are arranged and deformed, and slip or twining occurs between the atomic structures. Through deforming the arrangement of atoms in lattices, the coordination number of the central atom changes. Since the crystal distribution in the Fe structure is isotropic, the coordination changes regions are also randomly distributed during the simulation. It can be seen from Figure 5 that the number of 14-coordinate and 13-coordinate atoms gradually increased during the simulation. Due to the formation of defects, such as dislocations, some iron atoms form a 16position coordination. The specific structure of different coordination numbers is shown in Figure 6. We can deduce that with the stretching process, the coordination number of Fe in the entire system will change. Therefore, the Fe phase is not easily broken and shows high tensile strength. In addition, we discussed the relationship between strain and number of coordinates.  When the strain reaches 0.04 Å/Å, the stress reaches the elastic limit. As shown in Figure 7, the coordination number does not change significantly until the strain reached 0.04 Å/Å. In elastic behavior, when the material is under load, the flexibility between the metal bonds does not cause a significant change in the coordination number. At the macro level, the material will not suffer substantial damage during elastic deformation. After entering the yield stage, the coordination number decreased significantly. After the strain reached 0.08 Å/Å, the coordination number showed a small upward trend, which is mainly because of dislocation strengthening in the yield stage. After entering the yield stage, the material undergoes yield strengthening and the mechanical properties will gradually increase. This discipline is consistent with macroscopic phenomena. In the elastic phase, the coordination number maintains constant. In the yielding stage, the coordination number decreased. After breaking, the coordination partial recovery. Therefore, the stretching process is essentially breaking the coordination of atoms in materials. The similar analysis is taken for γ-Fe 2 O 3 phase and γ-FeOOH phase.
As shown in Figure 8, for the structure of γ-Fe 2 O 3 , the coordination number curve showed a significant decrease when the strain reached 0.08 Å/Å. Corresponding to the stress-strain curve, the material passes the elastic limit and reaches the yielding stage after the strain reaches 0.08 Å/Å. It can be observed from Figure 10 that the tensile failure of γ-Fe 2 O 3 starts form both ends and gradually moves toward the middle of the structure as the simulation time extends. The phenomenon is caused by periodic boundary conditions set during the simulation. During the stretching process, the Fe-O s bond (as shown in Figure 9) evolved from the initial stable 6 coordination to the 5, 4 coordination, and a small amount evolved to the 3 coordinates. Where the damage occurred, the  coordination is much lower than the average. As the damage progresses further, the distribution of atomic coordination changes gradually becomes disordered. Comparing to the Fe phase, in the γ-Fe 2 O 3 phase system, the coordination structure is only changed in a part of region. The external force can't be dispersed throughout the whole system. After the breaking, the coordination number will not be recovered. Therefore, the oxidation can increase the brittleness of Fe phase. For explaining the different stretching process of γ-FeOOH at difference direction, we further studied the structural evolutions of γ-FeOOH stretching. When these materials are stretched at the same tensile rate, the mechanical properties exhibited by γ-FeOOH are greatly different when stretched in different directions, as shown in Figures 11, 12. The tensile strength at x direction is stronger than that at z or y direction, which means the mechanical properties of FeOOH is stronger in this direction. Destruction of the bond of Fe-O s in the x direction results in the crack of material. The coordination number indicates that the coordination distribution of Fe-O s in γ-FeOOH is relatively uniform. When stretched along the y direction, material damage is primarily achieved by breaking the H bonds formed by the hydroxyl groups between the layers and the layers. This phenomenon proof that γ-FeOOH performed different stretching properties at certain directions. In the x direction, the evolution of coordination is like the Fe phase according to the Figure 13. With the stretching process, the coordination number is changed in whole system. The external force is dispersed throughout the whole system. But, in the z direction, the evolution of coordination is like the γ-Fe 2 O 3 phase. The coordination structure is only changed in a part of region. The external force is not dispersed throughout the whole system. Therefore, the brittleness of γ-FeOOH phase is increased in z direction.  Frontiers in Materials | www.frontiersin.org September 2020 | Volume 7 | Article 523126 8 CONCLUSION In this work, the reactive MD simulation was employed to studying the mechanical properties of steel and its' passivation film. The uniaxial stretching of different compounds (γ-FeOOH, γ-Fe 2 O 3 , and iron) was performed. We found that oxidation can reduce the strength of iron, but increased ductility under certain conditions. The evolutions of micro-scale structure in the  Frontiers in Materials | www.frontiersin.org September 2020 | Volume 7 | Article 523126 9 destruction process of each compounds were clarified. The detail analysis revealing the relationship of internal bonding and mechanical properties of steel. The detail conclusion is as follow: Firstly, the oxidation can reduce the tensile strength of steel. For the three compounds of γ-FeOOH, γ-Fe 2 O 3, and Fe, the order of tensile strength from high to low is Fe > γ-Fe 2 O 3 > γ-FeOOH. But, the ductility of γ-FeOOH under x direction is increased. Moreover, in the stretching process, the Fe and γ-Fe 2 O 3 phase resist the stretching by the metal bonds or covalent bonds among atoms. But the γ-FeOOH resist the stretching by the chemical bonds and intermolecular interaction at x and z direction, respectively. The force between the layers of γ-FeOOH is weak and more susceptible to damage. Finally, the tensile strength is origin from the coordination in the materials. There are two kinds of stretching process. One is in the Fe phase and x direction of γ-FeOOH phase, in which the coordination number is changed in whole system with stretching. The external force is dispersed throughout the whole system. Other is in the Fe 2 O 3 phase and z direction of γ-FeOOH phase. The coordination structure is only changed in a part of region. The external force is not dispersed throughout the whole system. The difference of stretching process result in large difference of tensile strength.
This investigation leads to possible new direction for studying the tensile strength of materials. In the future, it will be desirable to reveal the fracture process of steel with composite structure of Fe, γ-Fe 2 O 3 , and γ-FeOOH in which there are different environmental conditions. These works of the reactive MD simulations provide an excellent approach for observing and revealing the stretching process of solid materials. This strategy of evaluating materials tensile strength can supply valuable information in evaluating and improving the mechanical properties of reinforced concrete.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the manuscript/Supplementary Material.

AUTHOR CONTRIBUTIONS
ZL, XX, and MW performed this work and are the primary contributors in this work. MW primarily wrote the manuscript. PW, XH, YZ, and XW performed this work in details and assisted in making the picture. All authors of this manuscript have reviewed the manuscript and approved it for publication.

FUNDING
Financial support from National Natural science foundation of China under Grant 51678317, 51978352, 51908308; Natural FIGURE 13 | Partial structure of Fe with different coordination in γ-FeOOH. Green balls represent oxygen on the hydroxyl group, and red balls represent oxygen atoms connecting two iron atoms.