Multiscale numerical simulation of in-plane mechanical properties of two-dimensional monolayers

Many applications of two dimensional (2D) materials are often achieved through strain engineering, which is directly dependent on their in-plane mechanical characteristics. Therefore, understanding the in-plane mechanical characteristics of the 2D monolayers becomes imperative. Nevertheless, direct experimental measurements of in-plane mechanical properties of 2D monolayers face great difficulties due to the issues related to the availability of high-quality 2D materials and sophisticated facilities. As an alternative, numerical simulation has the potential to theoretically predict such properties. This review presents some recent progress in numerically exploring the in-plane mechanical properties of 2D materials, including first-principles density functional theory, force-field based classical molecular dynamics, and the finite-element method. The relevant case studies are provided to describe the applications of these methods along with their pros and cons. We hope that the multiscale simulation methods discussed in this review will inspire new ideas and boost further advances of the computational study on the in-plane mechanical properties of 2D materials.


Introduction
The discovery of graphene through the mechanical exfoliation of bulk graphite opens the eld of two-dimensional (2D) materials. 1 To date, many 2D materials including hexagonal boron-nitride (h-BN), phosphorene, MXene, layered double hydroxide (LDH), metal-organic frameworks (MOFs), covalent organic frameworks (COFs), and transition metal dichalcogenides (TMDs) have been experimentally and theoretically investigated.  The mechanical properties of 2D materials have become a hot research topic because mechanical properties are physical features when a material is acted upon by external forces during the large-scale practical applications. 30 2D materials can be deformed by in-plane stretching or by out-of-plane bending. Therefore, the mechanical characteristics of 2D structures include both in-plane and bending moduli. For 2D monolayers, the in-plane mechanical properties have attracted specic attention since they are directly relevant to the applications through strain engineering, as illustrated in Fig. 1. [31][32][33][34][35][36][37][38][39] As outlined, the 2D materials are promising candidates for stretchable transparent electrodes with high conductance and high exibility. 40 The resistance of the graphene electrode displayed no evident variation up to a bending radius of 2.3 mm corresponding to the uniaxial tensile strain of 6.5%. Meanwhile, the pre-strained substrates could enhance the bending limit up to $11%. 2D nanoelectromechanical systems (NEMS) resonators are strongly dependent on their mechanical properties, e.g. Young's modulus, mass density, and resonant frequency. 41 The presence of piezoelectricity characteristic coupled with the mechanical exibility of some 2D materials enables their promising applications in wearable power-generated nano--devices. The associated piezoresistive and piezoelectric effects under mechanical strain in 2D structures extend their applications to a diverse range in nano-industries. [42][43][44][45] The 2D hybrid organic-inorganic perovskites (HOIPs) solar cells possess power-conversion efficiencies surpassing 25%. 46,47 These 2D structures are promising applicants for the next generation of optoelectrical nano-devices, which bandgap can be tuned through strain engineering. 48,49 Jiao et al. conducted a systematic study to investigate the structural and mechanical properties of TMDs using DFT calculations. Their theoretical studies revealed that the 2D TcS 2 and TcSe 2 monolayers exhibit promise as potential candidates for light harvesting. 50 Furthermore, the analysis of the band alignment relative to the vacuum level showed that the TcSe 2 monolayer is potentially plausible for water splitting. Meanwhile, a 2% compressive strain can also make the TcS 2 monolayer suitable for photocatalysts. Considering the importance of strain engineering for 2D materials, many efforts have been devoted to exploring their in-plane mechanical characteristics. Table  1 lists the details of typical in-plane mechanical properties under the different external forces. The in-plane elastic moduli, Young's moduli, shear moduli, and Poisson's ratio are the parameters widely investigated in most studies.
There are different kinds of experimental technologies developed to measure the in-plane mechanical properties of 2D materials. These include (1) nanoindentation of suspended monolayer using atomic force microscopy (AFM); and (2) pressurized blister tests. Lee et al. rst used the AFM nanoindentation method to investigate the mechanical characteristics of graphene. 51 They reported that the forcedisplacement properties were interpreted within a framework of nonlinear elastic stress-strain response and yields an elastic stiffness of 340 N m À1 . Such quantities reveal graphene as comparatively the strongest and toughest materials ever investigated. Regardless, the experimental measurements require highly sophisticated facilities. Additionally, it is challenging to synthesize large high-quality 2D monolayers for the analysis. As such, numerical simulations to investigate the mechanical characteristics of 2D materials has become attractive, as evidenced by the blue area in Fig. 2.
Computational materials science has become a rapidly growing multidisciplinary eld. Striving to advance computing capabilities to comprehend, understand, and solve the complex problems of functional materials, including their mechanical properties. The mechanical properties are derivable from the linear stress-strain relationship or quadratic energy-strain relationship. Computational approaches with different size and time scales have been developed to explore the microscopic and macroscopic mechanical behaviours of 2D materials. The methods range from the sub-atomic [e.g. rst-principles density functional theory (DFT)], atomistic level [e.g. force eld-based classical molecular dynamics (MD)] and macroscopic levels [e.g. nite element analysis (FEA)] for the process simulation and engineering design, as illustrated in Fig. 3. In this review, the principles of the numerical simulation of the in-plane mechanical properties of 2D materials are presented. Some recent cases facilitate further discussion about the pros and cons of different methods.

DFT method
The relationship between the in-plane elastic constants and moduli of 2D monolayer materials can be described based on Hooke's law under the in-plane stress condition. 2 4 s xx s yy s xy Here C ij (i,j ¼ 1,2, and 6) is the in-plane stiffness tensor or inplane elastic constants, which are equal to the second partial derivative of energy (E s ), with respect to strain (3). For a rst-principles DFT method, the in-plane mechanical properties of 2D materials are calculated through the parabolic dependence of the energy on the elongation. [52][53][54][55][56][57] The energy (E s ) can be obtained from the DFT calculations with the different intervals of elongation. Using high-symmetry graphene as an example: the uniaxial strain 3 was applied along the x or y direction, which leads to 3 yy ¼ 0 or 3 xx ¼ 0, respectively. Under this condition, C 11 and C 22 were calculated from the coefficient of the quadratic term by tting the data of elastic strain energy [E s (3)] as a function of strain (3). Then, C 12 was calculated based on the equation while the equi-biaxial strain is applied. Finally, C 66 was calculated by tting the second-order polynomial.
Another method to calculate the in-plane mechanical properties is through the standard approach. In this method, the conventional mechanical properties are rst calculated. Aer that, the multiplication of the conventional stiffness tensors and monolayer thickness will equate to the in-plane stiffness tensors. Here, the thickness of the monolayer needs to be carefully justi-ed. Additionally, the vacuum used to separate the monolayers under the periodic boundary conditions also needs to be thick enough for the converged in-plane stiffness tensors to be retrieved.
The in-plane planar Young's (Y 2D ) and shear (G 2D ) moduli, together with Poisson's ratio (n 2D ) of 2D monolayer materials, can be derived from the in-plane planar elastic constants as:

MD method
Each atom for MD simulation is subjected to Newton's law of mechanics. The sequential steps to conduct MD simulation can be presented by the owchart illustrated in Fig. 4.

Stress-strain method.
Investigating the in-plane mechanical properties (e.g. Young's modulus) of 2D materials required the uniaxial tensile test to obtain the stress-strain relation. 58 The initial step for the stress calculation is to dene the internal pressure tensor P for N particles in volume V: where m i and v i are the mass and the velocity of the i th particle. 5 symbolizes an outer product and W is a virial tensor. The interaction between the atoms are accounted for by the force eld potential dened as: 59 In this equation: N is the total number of atoms, i and j are atomic indices, and u ij is the empirical potential. u ij is used to describe the interaction energy between atoms. Each energy term u ij depends only on a small number (N k ) of atoms located at positions r 1 ,r 2 ,.,r Nk , including atoms i and j. These atoms can be either in the supercell or in periodic images. Groups are assumingly chosen so that i is always in the supercell, which is denoted as 0. A periodic image of the group belongs to the periodic cell n. As such, the potential energy can be expressed as a sum of elements of the primary cell in an unambiguous manner: The partial force on a specic atom image can be calculated by deriving this expression with respect to the position r in of the periodic image of atom i in cell n: Fig. 4 The principle approach of molecular dynamics simulation.
According to the theory of Thompson et al., the virial tensor is: 60 The true mechanical stress tensor is then simply the opposite of the internal pressure tensor: The stress is proportional to the strain when the strain is small, which obeys the general Hooke's law. The Young's modulus can then be calculated because it is the slope of the strain-stress curve, as illustrated in Fig. 5. 61 2.2.2 Energy-strain method. Identical to the DFT method, the in-plane mechanical properties can be calculated based on the quadratic energy-strain relationship. In the MD simulation, the energy E s is calculated using the empirical force eld potential based on the atomic congurations of 2D monolayers. The local environment includes the inner state of a molecule, such as the bond length, bond angle, electrostatic forces and van der Waals forces. The total energy of a molecule system is the sum of its kinetic energy and potential energy. Additionally, the force eld is a function of potential energy which, can be expressed by the following equations: where E total refers to total energy of the system, E valance is bonding energy, E nonbond is non-bonding energy, E bond is bonding stretching term, E angle is angle bending term, E torsion is torsion angle term, E inversion is inverted angle term, E vdW is energy generated by van der Waals force, E electrostatic is energy generated by coulomb electrostatic force and E Hbond is hydrogen bond energy. Different types of energy terms are shown in Fig. 6.

FE method
The FE is a numerical method for solving problems, which can be described by partial differential equations or formulated as functional minimization. A domain of interest is represented as an assembly of nite elements. The FEA exhibits two main features. Firstly, the piece-wise approximation of physical elds on nite element provides good precision even with simple approximating functions. Another characterization of the approach is that the locality of the approximation leads to spare equation systems for a discretized problem. Thus, helping to solve problems with a massive number of nodal unknowns. The FE modelling approach was rst utilized in 2003, where the theory of classical structural mechanics was extended into the modelling of nanostructures. 62 Atoms are bonded together with a corresponding length and angle in a three-dimensional (3D) space. When subjected to loading, materials behave like space frame structures. Thus, the bonds between atoms are considered as connecting load-carrying generalized beam members. Meanwhile, the atoms act as joints of the members, as illustrated in Fig. 7. Structural mechanics analysis is conducted to determine the displacements, strains, and stresses of a structure under given loading conditions. The stiffness matrix approach is one of the elementary instances for solving linear elastic static problems. An additional application is to solve problems that deal with mechanical characteristics such as buckling, plasticity and dynamics. The stiffness matrix C consists of the following submatrices: where  The above elemental stiffness matrices indicate four parameters that are needed including the tensile resistance EA, the torsional stiffness GJ, and the exural rigidity EI x and EI y . All the while, identifying the length of the element, L. The EA, GJ and EI can be determined based on the energy equivalence since the deformation of a space frame results in the changes of strain energies. Based on the theory of classical structural mechanics, the strain energy of a uniform beam of length L subjected to pure axial force N can be presented as following: where DL is the axial stretching deformation, see Fig. 8a.
The strain energy of a uniform beam under pure bending moment M is Fig. 7 Simulation of a graphene sheet as a space-frame structure. Fig. 8 Pure tension, bending, and torsion of an element.
where a denotes the rotational angle at the ends of the beam, see Fig. 8b. The strain energy of a uniform beam under pure torsion T is where Db is the relative rotation between the ends of the beam, see Fig. 8c.
On the other hand, this potential energy between atoms is in the sum of contributions from bond stretch interaction U r , bond angle bending U q , dihedral angle torsion U f , improper (out-of-plane) torsion U u , and a non-bonded van der Waals interaction. 63 For covalent systems, the main contributions to the total steric energy come from the rst four terms, which have included four-body potentials. The harmonic approximation is adequate for describing the energy, considering the assumption of small deformation. 64 The dihedral angle torsion and the improper torsion are oen merged into a single equivalent term for convenience and simplicity's sake, i.e., where spring constants k r , k q , and k s are the bond stretching force constant, bond angle bending force constant and torsional resistance respectively. Meanwhile, the symbols Dr, Dq, Df represent the bond stretching increment, the bond angle change, and the angle change of bond twisting, respectively.
Here, both U r and U A represent the stretching energy, both U q and U M represent the bending energy, and both U s and U T represent the torsional energy. It can also be assumed that the rotation angle 2a is equivalent to the total change Dq of the bond angle, DL is equivalent to Dr, and Db is equivalent to Df. Hence, by comparing eqn (17)-(22), a direct relationship between the structural mechanics parameters EA, EI and GJ can then be deduced from the molecular mechanics parameters k r , k q and k s as following: Applying this concept followed by the solution procedure of stiffness matrix method for frame structures, the deformation and related elastic behaviour of 2D materials can be simulated.

Applications of numerical simulations
3.1 DFT studies 3.1.1 Graphene. Investigating the elastic properties (e.g. Young's and shear moduli and Poisson's ratio) of 2D materials reveals signicant behaviour of the systems and characterizes their structural stability. Various DFT-based studies with differing exchange-correlation functionals have identied the in-plane mechanical properties of pristine graphene. [65][66][67][68][69][70][71] Two different sets of units: (1) GPa and (2) N m À1 are oen used in the study of 2D monolayers. GPa is for the conventional mechanical properties, and N m À1 is for the 2D mechanical modulus, which is converted to the conventional unit through the division of monolayer's thickness. For example, the calculated 2D Young's moduli by Kudin et al. 72 is 345 N m À1 . Then, its conventional Young's modulus is 1.03 TPa by assuming the think of the graphene is 3.35 A, which is in good agreement with the measured values. We have used the DFT-D3 method to calculate the in-plane stiffness of graphene. The in-plane elastic constants C 11 , C 12 and C 66 are 348, 80 and 133 N m À1 , respectively. Based on eqn (5) and (6), the calculated in-plane 2D Young's modulus and Poisson's ratio is 329 N m À1 and 0.23, respectively. It Implies that the choice of exchange-correlation functional is not sensitive to the calculated 2D mechanical properties of graphene. Likely ascribing to the high in-plane stiffness of graphene, thus directly resulting in a hexagonal lattice and strong carbon-carbon covalent bonds.
3.1.2 TMDs. Some 2D materials can display several atomic structures that are either sufficiently close in energy to be experimentally observable or stabilized through substrates or doping. For such 2D congurations, the atomic arrangement can signicantly affect the physical properties of the system 73-91 (e.g. TMDs). TMDs are one family of well-studied 2D material due to their diverse physical properties in terms of their structural phases. 83,92 As an instance, molybdenum disulde (MoS 2 ) is reported to be most stable in its 1H semiconducting-like structural phase. However, it has also displayed to be metallic and unstable in its 1T structural phase. Additionally, it exhibits a small bandgap topological insulator in its 1T 0 phase. 74,75 The formation of their lateral heterostructures allows a new degree of exibility in engineering electronic and optoelectronic devices. Imani Yengejeh et al. carried out a comprehensive DFT calculation to study the impact of the structural phases on the mechanical properties of 1H, 1T 0 and 1H/1T 0 heterostructure phases of different TMD monolayers. Including MoS 2 , molybdenum diselenide (MoSe 2 ), tungsten disulde (WS 2 ), and tungsten diselenide (WSe 2 ). 83 Their results reveal that the impact of the structural phase is signicant (Fig. 9). The elastic constants of 1H-MoS 2 were C 11 ¼ 134, C 22 ¼ 134, and C 12 ¼ 29 N m À1 which was in good agreement with the ones obtained in the literature (where C 11 ¼ 130, C 22 ¼ 130, and C 12 ¼ 32 N m À1 , respectively). 93 The lateral heterostructures had a relatively weak mechanical strength for all the TMD monolayers. As a result, the signicant correlation between the mechanical properties of the TMD monolayers and their structural phases can be used to tune their stiffness of the materials.
3.1.3 MXenes. Recently discovered MXenes are a wellknown category of chemically active 2D materials oen synthesised through the surface functionalisation using terminal groups (e.g. -OH, -O, -Cl, -F, or their combinations). 30 The resulting surface functionalisation led to numerous technological applications, e.g. energy storage, 3,94-97 electromagnetic interference shielding, 98-102 composite materials, 103-107 catalysts 108-111 and sensors. [112][113][114][115] The strain-tunable electrochemical properties of MXenes enable them to be a propitious solution for exible and stretchable devices. 3,116,117 Despite over 3000 publications to date on MXenes since their discovery, only 4.9% of the investigations were conducted on the mechanical characterisation of the 2D materials. 30 On the other hand, a few experimental studies on the mechanical stiffness and strength of MXenes have been published, ascribing to the challenges of the experimental procedures. 118,119 Therefore, the necessity of conducting more computational studies regarding the mechanical characterisation of MXenes has been increasingly highlighted.
Luo et al. 120 performed a comprehensive DFT-based study to investigate the mechanical properties of the functionalized MXenes. They calculated the elastic constants, free energy, and work function of different MXenes to explore and comprehend the impact of the transition metal and surface functional groups on their characteristics. Apart from some O-termination group exceptions, most of their material models were observed to be dynamically stable. They nally indicated that the Young's moduli for the most investigated MXenes uctuated from 150 to 400 N m À1 and concluded that the functionalized MXenes usually exhibits stronger structural stability. Kazemi et al. recently conducted a theoretical investigation to study the mechanical properties of 2D titanium carbide MXene-based materials utilizing DFT calculations. 121 They discovered that in-plane elastic constants, Young's and shear moduli increase over the thickness of the MXene systems. It further suggests that the improved mechanical properties of the functionalized MXenes. The stronger interaction between the -O terminal group and the Ti n+1 C n leads to the stronger mechanical properties of -O functionalized MXenes. Consequently, Ti 4 C 3 O 2 is the strongest material among all their considered systems. More importantly, all Ti 4 C 3 T 2 MXenes have the higher in-plane Young's and shear moduli than graphene, which was previously thought to be the strongest known 2D material. Therefore, Ti 4 C 3 T 2 can be a promising alternative to graphene for applications that need 2D materials with considerably high stiffness and largely resistant to shape change. A similar role for the surface functionalization was also reported by Fu et al. 122 Amongst the ve functional groups that they considered, (-F, -Cl, -OH, -H, and -O) the oxygen functionalized Ti 3 C 2 O 2 possesses the largest adsorption energy, the highest in-plane planar elastic modulus, and the greatest enhancement of strength. Those ndings pave the avenue for tuning the mechanical properties of MXene-based materials database by engineering their composition.
3.1.4 Other 2D materials. Up to date, many 2D monolayers have been investigated using the DFT-based energy-strain relationship. Table 2 lists some recent studies for evaluating the mechanical properties of 2D monolayers with different exchange-correlation (XC) functionals.
For example, Aghdasi et al. performed DFT calculations to study the in-plane Young's modulus of monolayer phosphorene with and without the adsorption of different atoms, including Li, Mg, O, Al, Pt, Pd, and Mo. 126 Using the energy-strain method, they found that the in-plane Young's modulus of monolayer phosphorene is anisotropic. Its in-plane Young's modulus along the armchair and zigzag directions are 25.36 and 92.30 N m À1 , respectively. They also found that the impacts of different adsorbed atoms on the mechanical properties can either be positive or negative. Additionally, the inuence of the atomic adsorption on the in-plane Young's modulus of the armchair phosphorene is more remarkable than the zigzag phosphorene. Observing the adsorbed zigzag phosphorene reached the plastic region at smaller strains.

MD simulations
MD based simulation can be a promising approach for investigating the physical and mechanical properties of 2D materials using a relatively large size due to the low computational cost.
3.2.1 Graphene. Both stress-strain and energy-strain methods have been used to study the in-plane mechanical properties of graphene. For the stress-strain method, tensile tests were conducted to identify the Young's modulus, the ultimate tensile strength, and the yield strength. The uniaxial tensile loadings enabled them to obtain the elastic modulus, fracture stress, and fracture strain of graphene with respect to the orientation, sheet size and sheet temperature. Lebedeva et al. 127 evaluated the impact of the different potentials on the in-plane mechanical properties of graphene using the energystrain method. The Tersoff potential (REBO-1990, REBO-2000, REBO-2002 and AIREBO) and reactive bond-order potentials (LCBOP, PPBE-G, ReaxFF-CHO and ReaxFF-C2013) were used in their study. They uncovered that the impact of the potential was signicant. Among all the force elds, LCBOP provided results consistent with the experimental and DFT data. Thus, describing the in-plane deformations of graphene. However, the ReaxFF potentials greatly overestimate the Poisson's ratio. Additionally, the elastic response of all the potentials considered is non-linear already at elongations of 3%. This nonlinearity, however, is particularly striking for REBO-2000, AIR-EBO and REBO2002. It is oen responsible for very different results for the Young's modulus and Poisson's ratio across varying intervals of strains and computational approaches.
Javvaji et al. used the stress-strain method to investigate the inuence of external parameters on the yield properties of graphene. 128 They also studied the distribution of the stressstrain plots with varying domain sizes, initial crack lengths and lattice orientations. The yield values were observed to drop signicantly for small initial crack lengths less than or equal to 0.1 L. The yield values were also observed to decrease with the increase in domain size, as the larger specimens offer more spaces for dislocations to initiate. Evidence was given for the combined effect of the domain size, lattice orientation and crack length on the yield values of stress and strain. Gamboa et al. conducted a MD simulation to investigate the elastic properties of graphene sheets through implementing the uniaxial tensile tests. 129 The values of the Young's modulus and Poisson's ratio obtained in their calculation were 730 GPa and 0.39, respectively. Differing signicantly from the former estimations and much closer to experimental values. A proposed explanation was that the effect of atomic relaxation leads to a plausible accuracy. They nally observed an extended linear domain in the stress-strain curves, which is relevant to Young's modulus calculations at nite temperature.
Anastasi et al. conducted an intriguing investigation on the mechanical properties of pristine and nano-porous graphene using MD simulations. 130 They studied the fracture behaviour, and the temperature dependence of the 2D material models. Additionally, they investigated the random and uniformly distributed vacancy defects. They discovered that higher temperatures tend to signicantly decrease the fracture stress and strain almost linearly with temperature. The elastic modulus was, however, only affected at system temperatures higher than approximately 900 K. Thus, concluding that the fracture stress decreases substantially inversely proportional to the defect density. Their latter attempt was made for potential applications which require non-pristine sheets such as ltration membranes.

TMDs.
A MD simulation with the stress-strain method was performed by Ying et al. 131 to study the mechanical behaviour of synthesized ternary TMDs under two loading conditions: the uniaxial tension along the armchair and zigzag directions and the biaxial tension along both directions. They examined the impact of loading direction, the temperature on the Young's modulus and also the fracture behaviour of their investigated TMD nanosheets. The force interactions between atoms in their system were described by the Stillinger-Weber (SW) potential, which parameters were achieved by tting the SW potential to an experimentally obtained phonon spectrum. Based on their MD results, the in-plane Young's modulus of monolayer MoS 2x Te 2(1Àx) was slightly affected when x is in the range of 0-0.4. Yet, ostensibly increased when x is larger than 0.4, which demonstrated that the Young's modulus of MoTe 2 is insensitive to doping S atoms inside them. Furthermore, it shows that the Young's modulus of MoS 2 nanosheets is very sensitive to the doping Te atoms inside the nanosheets. Moreover, MoS 2x Te 2(1Àx) was found to possess a ductile fracture feature. Their Young's modulus and ultimate strength decreased when the temperature became higher, which was ascribed to the temperature-induced soening effect.
Zhao et al. conducted a MD study to evaluate the temperature-dependent mechanical properties of monolayer MoS 2 . 132 Their computational results stated that the mechanical characterizations of the material model vary with temperature. More precisely, the Young's modulus, maximum loading strain, and maximum load stress decrease with the rise in temperature from 4.2 K to 500 K. They also discussed that the tendency of the maximum loading strain with the different temperature is opposite to that of metal materials which are caused by the short-range SW potentials among the TMD atoms. Despite holding the high computational efficiency, the mathematical form of the SW potential may not be suitable for the atomicthick planar structures, such as graphene and b-BN, as such potentials are less probably to resist the bending motion of these real planar crystals. [133][134][135] Jiang performed a MD simulation to study the mist straininduced buckling of the TMDs lateral heterostructures, such as MoS 2 -WSe 2 and MoS 2 -MoTe 2 . 136 He also utilized the SW potential to describe the interatomic interactions for each TMD and the coupling between different TMDs. Computational results indicated that mist strain presumably causes buckling of the TMD with larger lattice constants in the lateral heterostructure, due to the atomic-thick nature of TMDs. These results raise many mechanical challenges for the structural stability of the low-dimensional materials.
3.2.3 Other 2D materials. Due to the limited availability of force eld potentials, the studies on other 2D materials are few. Table 3 summarizes some recent studies for evaluating the mechanical properties of 2D materials using MD methods.
One recent MD simulation example conducted by Yang et al. investigated the temperature-dependent stress-strain relations of monolayer black phosphorus while considering the impact of temperature changes. 142 The Young's and shear moduli under uniaxial tension at different temperatures were investigated. The values of the Young's and shear moduli ranged from 24 to 19.2 GPa and 25.4 to 20 GPa, respectively. Observations showed that for a given temperature, the Young's moduli along the zigzag direction are more than four times higher than those along the armchair direction, which matches the DFT conclusion. 126 Their predicted strength and moduli were in good agreement with the available experimental data.

FEM simulations
Due to its signicant advantages, the FEM approach has become increasingly popular and a promising alternative for evaluating the macroscopic mechanical properties of lowdimensional materials for a wide range of emerging applications where high conductivity, mechanical strength, high surface area, and high yield strain are required. 143,144 FEM's features enable researchers to work on a larger scale of the material models, which consists of numerous atoms in their structural congurations.
3.3.1 Graphene. Scarpa et al. proposed a FE model and an approach based on cellular material mechanics theory to study the in-plane linear elastic properties of the monolayer graphene sheets. In their material model, 145 the C-C bonds were represented by equivalent mechanical beams having full stretching, hinging, bending and deep shear beam deformation mechanisms. Finally, using equivalent mechanical C-C bond characteristics, the inplane Young's, shear moduli and the Poisson's ratio were derived.
FEA is also strongly benecial for conducting the simulations to explore the mechanical properties of defective lowdimensional systems. Recently, Nakanishi et al. investigated the macroscopic mechanical properties of hollow-wall graphene gyroids by performing indentation tests with suitable interpretation by FE simulations. 146 They carried out two sets of investigations. Firstly, they used periodic cell calculations to explore the mechanical properties of their systems in terms of the relative density and cell wall characterization of the lattice. Then, the indentation simulations of a continuum with the effective properties of the gyroid were performed. They nally concluded that hollow-wall graphene gyroids combine sizedependent mechanical and electrical properties with a topology of high structural efficiency.
3.3.2 TMDs. Recently, Li et al. combined the AFM and FEM to investigate the in-plane Young's modulus of 1H/2H MoS 2 monolayer and bilayer. 147 The bimodal AFM is rst employed to gure out the effective spring constant between the microscope tip and sample. Aer that, the FEM is developed to quantitatively ascertain the effect of substrate stiffness on deformation. Using this combined method, they calculated the conventional in-plane Young's modulus of monolayer MoS 2 of 265 GPa aer removing the impact of the substrate. When assuming a 1H MoS 2 monolayer thickness of 6.15Å, the corresponding 2D inplane Young's modulus is 163 N m À1 and comparably close to the DFT value of 165 N m À1 . 148 Therefore, this combined method provides a convenient approach to calculate the inplane Young's modulus of 2D materials on a substrate.
3.3.3 Other 2D materials. Up to date, the FE studies on the in-plane mechanical properties of other 2D materials are scarce due to the restriction of the available spring constants. Table 4 summarizes the most recent studies for evaluating the mechanical properties of 2D materials using the FEM.

Discussion
Overall, numerical modelling approaches with different size and time scales have been developed for evaluating the in-plane mechanical characteristics of 2D monolayers.
The in-plane elastic properties of some of the most studied 2D materials are presented in Table 5. 51,83,121,[147][148][149][153][154][155][156][157] Basically, in-plane mechanical properties (such as Young's modulus, shear modulus, hardness, and Poisson's ratio) can be calculated through the linear stress-strain or quadratic energystrain relationships. DFT methods are independent of any empirical parameters, which provide the highest accuracy and exibility, as evidenced by the data listed in Table 5. The DFT inplane Young's modulus and Poisson's ratio of graphene and MoS 2 are almost identical to the measured values. There is a considerably large difference in the MXene system. Because materials used for experimental measurements have hybrid  139 Boron nitride-carbon nanosheet Mechanical properties Investigating the tensile characteristics of single layer BN-C nanosheets. It was shown that the BN-C nanosheet with parallel arrangement exhibits slightly improved mechanical resistance than the BN-C nanosheet with series arrangement Hou and Yang 58 Graphene oxides Tensile properties Carrying out tensile test to investigate the mechanical properties of graphene oxide sheets Javvaji et al. 128 Graphene Yield stress and strain Investigating the combined effect of domain size, crack length, and lattice orientation on the mechanical properties of graphene Anastasi et al. 130 Graphene Mechanical properties and fracture behaviour Carrying out uniaxial tensile loading to investigate the mechanical properties of pristine and nanoporous graphene. An increase in system temperature results in a signicant reduction in the fracture stress and strain Siriwardane et al. 140 Sulfur-functionalized MXenes Structural, stability, and ion dynamic properties Performing computational calculations at 400 K and bondlength analysis to study the physical properties of functionalized MXenes Cellini et al. 141 Diamond boron-nitride Mechanical properties and pressure induce phase transition Investigating the mechanical properties and pressure-induced formation of 2D diamond boron nitride. The results show a metastable local rearrangement of the h-BN atoms into diamond crystal clusters when increasing the indentation pressure terminal groups, this consequently leads to the theoryexperiment discrepancy. The simulation quality of MD simulation is greatly affected by the force elds used in the study, as evidenced by Lebedeva et al. on graphene. 127 Table 5 lists one of the best MD results, which are also close to the experimental data. These results are caused by the potentials for graphene and TMDs being purposely optimized. Up to date, graphene is the most studied system via the FEA method. However, its calculated in-plane mechanical properties are far from the experiments. For example, the calculated Poisson's ratio of graphene is about four times higher than the DFT and experimental values.
In addition, the rst principles DFT method does not require any pre-set parameters. Its simulations can study chemical reactions involving charge transfer, bond formation and cleavage with high accuracy. Thus, DFT-based simulations provide greater exibility for calculating the in-plane mechanical properties of 2D materials. As a comparison, only limited force elds are available for most of the 2D materials beyond graphene. The parametrization of the force eld either derived from the experimental data or rst-principles DFT results is a very time-consuming task. As such, MD-based methods cannot be widely used. The FEA simulation face the same problem because its quality is based on the spring constants, On the other hand, the force-eld-based classical MD methods use empirical interatomic potentials to describe their interaction energies. Thus, signicantly reducing the computational time using low cost. Consequently, the MD method can assess the mechanical properties of large 2D materials. Additionally, the impact of temperature on the mechanical properties can also be evaluated. The FEM has the lowest computational cost. Most of the FEM computations can be conducted using desktop computers. Consequently, it can be used to investigate some macroscopic mechanical properties.
To obtain the holistic picture of in-plane mechanical properties of 2D materials, a combination of different computational methods is a plausible approach by taking advantage of each one's strong features. So far, there have been some attempts to combine the computational methods and propose a rened technique to investigate the mechanical properties of the nanostructures. Imani Yengejeh et al. proposed a rened FEA to evaluate the vibrational properties of the defective CNTs and their modications. 158 They implemented the well-established FEA for the perfect material models. For appraising the properties of the defective structure, the accurate DFT structural relaxation was used as FEA cannot consider the rearrangement of the atoms aer the introduction of the defect (e.g. vacancy defect). Their obtained results were in good agreement with the experimental ndings, which suggested that such a combined multiscale modelling approach is feasible for investigating the mechanical characteristics of the 2D materials.

Summary and outlook
Over the past few decades, the family of low-dimensional materials has rapidly grown beyond graphene. Many novel 2D materials including, TMDs, MXene and phosphorene, have become signicantly important. The remarkable properties of 2D materials pave a bright future for basic research and largescale industrial applications. Due to several limitations of the experimental measurements, numerical simulation becomes a promising alternative to evaluate the in-plane mechanical properties of these 2D monolayers.
This review described some widely used computational approaches for investigating the in-plane mechanical properties of 2D material. Their features are discussed using some recent study results as examples. DFT-based computation is the most widely used approach to understand the in-plane mechanical properties of 2D materials since their results are reliable. All 2D monolayers can theoretically be investigated using the DFT method aer their atomic congurations are correctly built. The DFT results can be used to guide their design and screen for 2D materials at the atomic level to achieve desired mechanical properties. A major disadvantage of the rst-principles calculations is that only a few thousand atoms can be simulated. Consequently, seriously limiting the direct comparisons to experimental studies. In experiments, the measurements typically occur on the micrometre length scale. For such larger sizes, classical MD and FEM simulations are desirable. The minuscule computational costs of atomistic simulations using empirical force eld potentials and spring constants enable the mechanical properties of some 2D systems to be predicted. Basically, MD simulation is one of the molecular simulation techniques which refer to a set of approaches a conducting computational experiment on molecule models. This simulation method can be classied into two major categories: microscale dened in the range of 0.1-10 nm and mesoscale dened from 10 nm to 100 nm. However, the key ingredient for MD and FEM simulations is an empirical interatomic potential and spring constants, which determines the availability and accuracy of MD and FEM simulations. Since different computational approaches exhibit their advantages and disadvantages, the most appropriate computational technique can be selected in terms of the type of mechanical properties of the 2D material, the required accuracy of the results, and the features of the applied method.
Finally, the development of novel technologies becomes essential to generate a comprehensive understanding of the inplane mechanical properties of the 2D materials for their applications based on strain engineering. Several promising methods are discussed below: (1) A combining approach may be an excellent solution to describe a wide range of 2D materials properties as it would enable the incorporation of positive features from each method. Table 5 In plane Young's modulus (Y 2D ), shear modulus (G 2D ) and Poisson's ratio (n 2D ) of 2D materials obtained by different computational approaches and experimental measurements For example, the machine learning method has been used to develop the force elds from the DFT results, which can greatly expand the application area of the MD and FEM simulations. [159][160][161][162] (2) The combination of the experimental measurements and the numerical simulations also holds great potentials. As suggested by the recent study, the experimental measurements can be used to provide accurate spring constants, which can enable the reliable in-plane Young's modules using the low-cost FEA approach. 163,164 (3) Machine learning technology can also be used to accelerate the advancement of this area. The DFT-based computational data can be used to build the database. The ensemble learning system can be built using a diverse set of base regressors/classiers. Construction is done through both traditional learning algorithms (such as support vector machine, decision tree, kernel ridge regression, Gaussian mixture regression) and deep learning algorithms to identify the relationship between the atomic ngerprints and their inplane mechanical properties. 165,166

Conflicts of interest
There are no conicts to declare.