Prediction and simulation of mechanical properties of borophene-reinforced epoxy nanocomposites using molecular dynamics and FEA

2.1 Construction and energy minimization

Using DGEBA and DETA molecules, an amorphous cell with a density of 1.20 g·cm⁻³ at room temperature and periodic boundary conditions was created. The system was subjected to 1,000 steps of energy reduction to attain the closest local minimum. A 2-ns isothermal–isobaric (NPT) dynamics with a time step of 1 fs at 298 K and 1 atm was conducted to accomplish the global minimum of potential energy. The Bredesen, Nosé–Hoover thermostat [24,25,26] and the Andersen barostat [27] were employed during the simulation to control the system’s temperature and pressure, respectively. In order to assess if the system had attained equilibrium, the temporal evolution of the energy and temperature of the system was tracked. Covalent bond formation: After the structure had reached equilibrium, the distances between reactive atom pairs were measured. Covalent bonds were then produced between the reactive atom pairs that were deemed to be suitable (those within crosslinking cut-off distance). After that, a sophisticated minimization task was performed using 1,000 iterations to reduce the load exerted on the system. This was done to give the molecules enough kinetic energy and increase the likelihood that they would fall within the reaction cut-off distance. The remaining reactive atoms’ distances were once again measured, and if possible, new covalent connections between them were made. After each stage of bond formation, the degree of cross-linking can be assessed, and step 3, which consisted of a cyclic combination of bonding and minimization dynamics, was repeated. As a result, a convergence of 65% was obtained. The boxes in the diagram stand in for the simulation boxes in the prior and subsequent crosslinking stages.

2.2 Force field and simulation procedure

A force field has been utilized, which is incorporated into “Materials Studio” (MS) [28]. Based on the results of COMPASS II and “force field techniques,” valence parameters were generated and validated using condensed phase properties in addition to empirical data obtained for molecules in isolation. A Coulombic function for electrostatic interactions, a 6–9 Lennard–Jones potential for Van der Waals interactions, and three categories of energy terms – bonded energy terms, cross-terms, and non-bonded energy terms – make up the COMPASS II force field. These terms are listed in equation (1).

where E b denotes covalent bond-stretching energy, E θ denotes bond-angle bending energy, E φ denotes torsion-angle rotation energy, E χ denotes out-of-plane energy, E q denotes electrostatic energy, E vdW denotes van der Waals energy, and E b , b x , E b , θ , E b , φ , E θ , φ , E θ , θ x , E θ , θ x , φ denote cross terms representing the energy due to the interaction among bonds.

When creating an amorphous cell, the goal density was kept constant at 1.20 g·cm⁻³, which is the same as the actual density of the epoxy DGEBA. A total of 1,040 molecules of DGEBA/DETA have been employed in tightly packed amorphous cells. The structure has been stabilized using geometric optimization and energy minimization. NPT ensembles have been used in the Forcite module to study dynamic stability. An increased external pressure is used for both dynamics and annealing modules in order to enhance the density and homogeneity of the material and hence increase its mechanical strength. In order to promote the rearrangement of the material’s microstructure and enhance its mechanical characteristics, annealing under high external pressure can assist the system maintain its equilibrium by reducing the mobility of the constituent particles. However, it should be noted that the application of pressure may also cause modifications in the material’s physical and chemical characteristics that might influence its overall stability, depending on the particulars of the system and the annealing circumstances.

In order to maintain volume and density, annealing has also been done using the NVT ensemble. Energy reduction and the dynamic run have both been completed. The control parameters for the Forcite module are depicted in Table 1, respectively.

The simulation was run at room temperature, or 25°C (298 K), under 0.3 GPa of external pressure. The generated amorphous cell’s anticipated density, computed employing the control parameters listed in Table 1 and shown in Figure 2, is 1.209 g·cm⁻³, which is quite close to the epoxy DGEBA’s actual density at room temperature. Figure 3a and b shows the temperature stability of the developed model over time. The temperature of the model first rose until it reached room temperature for NPT, as seen in Figure 3a and b. As the simulation continues, the temperature profile remains constant, indicating that the model is dynamically stable. Furthermore, the Forcite Dynamics energies trajectory has been validated by the stability of the developed model. The model was cooled down at the same ensemble assumptions.

Figure 2

Simulated data for density (g·cm⁻³).

Figure 3

(a) Simulation data for temperature NVT. (b) Simulation data for temperature NPT.

As the simulation goes on, it can be seen from Figure 4 that all energies – potential energy, kinetic energy, non-bond energy, and total energy – of the epoxy model initially rise and oscillate with the rise in temperature. The oscillation of the graph is observed up to 5 ps and then the model stabilizes and the energy remains constant throughout 10 ps. In Figure 4, as the simulation advances, it is observed that the total energies remain constant, and a decrease in kinetic energy leads to a more stable model, which also supports the model’s thermal stability.

Figure 4

Forcite energy diagram for the simulated amorphous cell of DGEBA.

For clarity, in Figure 5, the Ball-Stick styles are used to represent DETA molecules in their cross-linked structure rather than their non-cross-linked counterpart, and vice versa.

Figure 5

Amorphous model of epoxy resin with DETA.

2.3 Modelling of borophene β₁₂ layered sheets

The library of MS has “BN” (boron) crystal structure. By replacing the C atoms in the unit cell by B atoms and by another geometry optimization, the unit cell is optimized. The cleaving process along certain planes is very important. A crystal was built using a vacuum slab tab. The final step is to optimize it again, i.e., geometry optimization by a best-fit force field. A single layer of boron known as borophene was built as shown in Figure 6.

Figure 6

Crystal structures of borophene sheets of 2D boron materials. Top views of β₁₂.

η-LB layers were built by super celling a single unit cell of β₁₂, a form of borophene as shown in Figure 6. The layers and geometry will be reconstructed, and the geometry will be optimized with the Dmol3+ module. A supercell form of one-layer borophene is shown in Figure 7.

Figure 7

Supercell structure of a borophene sheet of 2D boron materials.

For β₁₂ orthorhombic borophene, the young modulus and the Poisson’s ratios across the x and y directions are computed in terms of four nonzero elastic stiffness constants. Using the control values indicated in Table 1, similar simulation techniques have been used. The created model was further processed for predicting mechanical properties once its stability was confirmed. The simulations were performed with temperatures fixed to 300 K in accordance with a set of strain rates from 0.04. The system was equilibrated utilizing the Nosé–Hoover thermostat (NPT) ensemble (fixed number of atoms, pressure, and temperature) for 5 ps. Table 2 exhibits the computed elastic stiffness constants, Young modulus, and Poisson’s ratios for η-LB.

2.4 Modelling of borophene/epoxy unit cell system

In this simulation, the polymer matrix is composed of DGEBA that has been cross-linked with DETA. As seen in Figure 8, the voids in the composites were filled with single- and multilayer borophene sheets. Chemical functionalization might be an excellent solution if you want to produce good graphene dispersion, while also producing a solid adhesive interface between the nearby polymer chains and graphene. Covalent bonds between functional groups, which provide chemical functionalization, are the foundation of chemistry. Before borophene may be employed as a successful reinforcement in polymer composites, there are many obstacles to be addressed. Among these difficulties is the dispersion of borophene, which is regarded as one of chemistry’s trickiest puzzles. Different functional groups are considered when talking about functionalized graphene; the amine group is highly reactive and can be directly incorporated into the epoxy resin [29]. In addition, various functional groups are examined in relation to functionalized graphene. The issue of borophene dispersion and aggregation in the epoxy matrix might be solved by using this method. In the epoxy matrix depicted in Figure 8(a), a single-layer borophene sheet containing 185 atoms has been randomly reinforced with 1,040 epoxy molecules. The cell dimension is 25.9 Å × 25.9 Å × 25.9 Å. In the current study, reinforcements have been carried out at various borophene weight percentages, including 5, 10, 15, and 20%. Figure 8(b) shows 1–4 layers stacked in between the epoxy, with an average borophene layer spacing seen to be approximately 2.5 Å. The cell dimensions changed with the change in layer stacking.

Figure 8

(a) Borophene sheet mixed in DGEBA/DETA matrix in a unit cell. (b) η-Layers of borophene sheet stacked in a unit cell.

2.5 Finite element analysis (FEA) of borophene/epoxy layer composite

ABAQUS 6.13 (FEA) software was used to analyse the layer composite in order to determine whether it has the rigidity and strength required for vessel building. In order to stack the laminate layers, the composite layers are arranged in various thicknesses and ply orientations. The composite vessel’s flexible design allows it to be adapted to a variety of working conditions by altering the ply stacking patterns and vessel shape characteristics.

At a given weight percent layer addition of borophene, the curve of the borophene-reinforced epoxy nanocomposite is nonlinear as shown in Figure 9(a–c). This may be due to the incorrect interphase or poor dispersion of borophene in a significant quantity.

Figure 9

(a) FEA analysis model. (b) Ply layup arrangement. (c) Meshed model.

FEA was performed on a composite sheet made of 2-LB β₁₂ borophene with a shell planar feature that was 110 mm × 50 mm in size. Material is assigned in the mechanical behaviour segment considering its mechanical property values from MD simulation. The values are assigned under the elastic lamina type. Under the material property tab, conventional shell element-type composite layup is assigned with a three-ply count and part global orientation with automated calculated section symmetry using the lamina property to simulate. The thickness was set to 1 mm, with integration points at 0°, 45°, 90°. For proper boundary conditions, an interaction step with two reference points was defined . A uniform distributed load of magnitude 100 N was defined for the model. The mesh size was kept global for proper and fast simulation, and the model was simulated with proper prior checks. The values assigned under the elastic lamina type are listed in Table 3.

3.1 MD simulation

On two amorphous cell models, MD simulation has been used to determine their properties. One model is made up of cured DGEBA/DETA with a borophene sheet layer mixed in percentage proportion. Another model is formed as a result of the incorporation of η-LB β₁₂ borophene sheet into the matrix of DGEBA. Specifically, the simulation was carried out in accordance with the considerations indicated in Tables 4 and 5. Using the Forcite dynamic simulation with continuous strain reduction, we were able to derive the mechanical characteristics. In the Forcite module, the simulations using the COMPASS II force field have been carried out using the data. At room temperature, the simulations have been carried out with the same parameter control as shown in Table 1 in the mechanical property module. From Tables 4 and 5, it is clear that the stiffness value of the epoxy increases significantly with escalating amounts of homogenous mixed η-LB β₁₂ reinforcement as the amount of reinforcement is increased. About three to four η-LB layers of β₁₂ borophene are required to get superior properties. When comparing 1-LB to other 2D monosheets, such as phosphorene, silicene, and MoS₂ monolayers [30–32], it is clear that the borophene layer exhibits much better mechanical attributes when considering the entire criterion. In the graph shown in Figure 10, the relationship between Young’s modulus and strain rate at various weight percentages of borophene is mapped. There is a significant rise in Young’s modulus with reinforcement at the same strain rate, which is comparable with the prior research findings.

Figure 10

Young’s modulus (GPa) versus strain rate for stacking up in layer form mixing of η-LB.

Figure 10 depicts a plot of Young’s modulus versus strain rate at various borophene weight percentages and without reinforcement. As shown in Figure 11, borophene reinforced epoxy nanocomposites show nonlinear behavior at a weight percent borophene layer addition. This may be caused by an incorrect interphase or poor dispersion of a significant quantity of borophene in the epoxy matrix. At 0.001 strain, the Young’s modulus of every borophene-reinforced epoxy structure (5, 10, 15, 20% weight percent reinforcement) at the same strain rate shows a significant rise in Young’s modulus, but it has also been reported that the Young’s modulus decreases with an increase in strain rate.

Figure 11

The Young’s modulus (GPa) versus strain rate varies as a function of the stacking up in layers form mixing of η-LB.

The modulus of elasticity of borophene epoxy is shown in Figure 11, reinforced with layer reinforcement of borophene at various strain rates. There is a considerable enhancement in the value of elasticity of the epoxy with a raising reinforcement value. The remarkable properties are obtained in 3-LB and 4-LB borophene. In addition, carbon nanotube (CNT)-reinforced epoxy nanocomposites and the findings of graphene-reinforced epoxy nanocomposites are compared [34]. When compared to CNT-reinforced epoxy nanocomposites, it has been discovered that graphene-reinforced epoxy nanocomposites display superior mechanical characteristics.

3.2 FEA simulation results

The microscale stress–strain curves that suggested fracture processes at the molecular level were successfully given by MS. The discrepancy between the ideal and real strength of materials should be carefully studied since fracture strength in the atomistic simulation can be understood as a type of ideal strength of materials. However, the MS analytical results, which undoubtedly discriminate between a matrix failure and an interface failure, will be extremely helpful in the material selection for the polymer composite design. In the simulation, the “encastre boundary condition” limits the structure of all “active degrees of freedom” within the defined area. This restriction will be applicable to all “nodes” within the defined area once the part is “meshed” and the task or assignment has been assigned.

The deformation and stress distribution at a given load are depicted in Figures 12(a and b) and 13(a and b). In this work, a unidirectional composite’s transverse failure properties are modelled. A transverse fracture, often known as a first-ply failure, causes significant harm and may result in deadly delamination. The composite has a significant capacity for plastic deformation, which can safely limit breakage failure under high pressures. The maximum Mises stress on the layer’s rising internal pressure is depicted up until it ruptures. This is owing to the composite ply’s enormous fibre failure at burst pressure, which causes high-stress concentrations and stress redistributions. In Figures 14 and 15, the stress–strain curve and load–displacement curve give an obvious result that agrees well with Yang et al.’s [35].

Figure 12

(a) Composite plate deformed under load condition. (b) Composite plate contours deformed under load condition.

Figure 13

(a) Composite plate material plot under deformed shape load condition. (b) Composite plate material arrangement under load condition.

Figure 14

Stress vs strain plot under load condition.

Figure 15

Force vs displacement plot under load conditions.

The findings of prior research investigated the mechanical properties of CNT-reinforced epoxy nanocomposites using MD simulations [34]. Such simulations can be used to analyse the behaviour of materials at the atomic and molecular scales and can help researchers understand how different factors, such as the concentration and alignment of CNTs, affect the mechanical properties of the composite [34]. The prior study has discussed the use of MD simulations and FEA to study the mechanical behaviour of oxygen-functionalized graphene/polymer nanocomposites [35]. In this case, the simulations would be used to study the behaviour of the graphene/polymer nanocomposites at the atomic scale, taking into account factors such as the concentration and alignment of the graphene particles and the implications of oxygen functionalization on the mechanical properties of the composite. FEA is a numerical method used to study the behaviour of systems with complex geometry. It is used to predict how a material will respond to specific loading conditions and can be used to study the mechanical behaviour of the oxygen-functionalized graphene/polymer nanocomposites by simulating the load applied to the material and predicting the resulting deformation and stress distribution. The outcomes have discussed the results of these simulations, such as the influence of oxygen functionalization on the mechanical properties of the composite and the relation between the concentration and alignment of the graphene particles and the mechanical behaviour of the nanocomposite [35]. The findings from the previous study reported vibration properties of foam plates that have been sandwiched and integrated with composite materials [36–38]. The study likely looked at how the addition of composite materials affected the damping and frequency of vibrations in the foam plates. They would have measured the changes in vibration properties at different frequencies and observed how the composite materials affected the overall damping of the foam plates [36,39,40]. There are studies focused on developing a molecular mechanics model to predict the buckling behaviour of multi-walled carbon nanotubes (MWCNTs) under axial compression [37,41,42]. The study likely used computational simulations to model the behaviour of MWCNTs under various levels of compression and compare the results to experimental data [42,43,44]. The goal of the study was likely to create a model that can accurately predict the buckling behaviour of MWCNTs under compression, which can be useful for understanding the mechanical properties of these materials and designing structures that utilize them [43–45]. The results of the study may have shown the developed model to be in good agreement with the experimental data and has the ability to predict the buckling behaviour of MWCNTs under axial compression in a more accurate way [37,46,47].

The study used computational simulations or mathematical analysis to investigate how the shape and material properties of the sheets affect their buckling behaviour under compression [47–49]. The goal of the study was likely to better understand the buckling behaviour of polygonal thin sheets and to provide design guidelines for structures that use these materials [50–52]. The results of the study may have shown the relationship between the shape of the polygon and the critical buckling load of the sheet, as well as the effect of different material properties on the buckling behaviour [53–55]. The study may have also presented a method to predict the buckling behaviour of polygonal thin sheets with different shapes, sizes, and materials [45,56,57]. The aim of the study was likely to understand how the addition of MWCNTs affects the elastic properties of high density polyethylene (HDPE) and how these properties are influenced by the MWCNTs’ loading and distribution in the HDPE matrix [58,59]. The results of the study may have shown that the addition of MWCNTs increases the elastic modulus, and the strength of HDPE may also have shown the effect of MWCNTs’ loading and distribution on the elastic properties of the composites. These results may be useful for understanding the mechanical properties of HDPE–MWCNT nanocomposites and designing structures that utilize these materials [60–62]. The goal of the study was likely to understand how thermal expansion and composite material properties affect the structural stability of carbon foam beams under axial compression and provide design guidelines for these types of structures [41,46,47]. The results of the study may have shown how thermal expansion and material properties affect the critical buckling and bending loads of the sandwich beams, as well as the deformation, and failure modes of the beams under compression and at elevated temperatures [41,48,49].