Nonlinear Static Bending and Forced Vibrations of Single-Layer MoS2 with Thermal Stress

Single-layer molybdenum disulfide (MoS2) has been a research focus in recent years owing to its extensive potential applications. However, how to model the mechanical properties of MoS2 is an open question. In this study, we investigate the nonlinear static bending and forced vibrations of MoS2, subjected to boundary axial and thermal stresses using modified plate theory with independent in-plane and out-of-plane stiffnesses. First, two nonlinear ordinary differential equations are obtained using the Galerkin method to represent the nonlinear vibrations of the first two symmetrical modes. Second, we analyze nonlinear static bending by neglecting the inertial and damping terms of the two equations. Finally, we explore nonlinear forced vibrations using the method of multiple scales for the first- and third-order modes, and their 1:3 internal resonance. The main results are as follows: (1) The thermal stress and the axial compressive stress reduce the MoS2 stiffness significantly. (2) The bifurcation points of the load at the low-frequency primary resonance are much smaller than those at high frequency under single-mode vibrations. (3) Temperature has a more remarkable influence on the higher-order mode than the lower-order mode under the 1:3 internal resonance.


Introduction
Since monolayer graphene was first mechanically exfoliated from graphite in 2004 [1], its excellent physical, chemical, and mechanical properties have attracted extensive attention [2][3][4][5][6][7].At the same time, the graphene-like two-dimensional (2D) transition metal dichalcogenides (TMDCs) have attracted widespread attention due to their single-layer characteristics and their excellent mechanical properties similar to those of graphene [8][9][10][11][12][13].Molybdenum disulfide (MoS 2 ) is a typical TMDC material, it can be obtained using mechanical stripping, a chemical approach, CVD synthesis, and other methods [14][15][16].There are significant differences in the size, quality, and yield of monolayer molybdenum disulfide prepared using different methods.MoS 2 not only overcomes the zero-band-gap drawback of graphene but also retains its numerous advantages.This makes it suitable for a broad range of potential applications [15,17,18].Thus far, research on MoS 2 has focused on its electrical, thermal, and friction properties [19], whereas its mechanical properties have rarely been investigated.TMDCs have been used as high-quality nanoresonators [20,21].Because the band structure of monolayer MoS 2 can be changed according to the mechanical strain, new nanomechanical devices can be designed by applying mechanical deformation.For example, Andres [22] fabricated a single-layered mechanical resonator using MoS 2 and demonstrated nonlinear behavior under room temperature and vacuum conditions.However, how to model the mechanical properties of 2D nanomaterials, in which the materials have a monolayer structure with single or multiple atoms, remains an open question.Because monolayer MoS 2 can resist bending deformations, the macroscopic Föppl-von where D = Eh 3 /12 1 − ν 2 is the bending stiffness; D 1 = Eh/ 1 − ν 2 is the extensional stiffness; and E, ν, and h are the elastic modulus, Poisson's ratio, and thickness of the plate, respectively.This indicates that the bending and extensional stiffnesses are in D/D 1 = h 2 /12 for classical plate theory, whereas the two stiffnesses in 2D monolayer nanomaterials are independent, namely D/D 1 ̸ = h 2 /12 [23].Therefore, one cannot obtain the out-of-plane bending and torsional stiffnesses using the in-plane mechanical parameters and the thickness of the 2D monolayer materials.This is called the Yakobson paradox [23], whereas some authors believe that there is no paradox because much of the literature fails to distinguish h and E from the effective thickness and the effective elastic modulus [24].Similarly, for a single-layer MoS 2 , one cannot directly obtain the out-of-plane bending and the torsional stiffness through its thickness.Based on the bond orbital theory of covalent bonds, Huang [25] obtained a continuous mechanical theory of monolayer graphene to explain the Yakobson paradox physically.This theory clarifies the physical mechanism of graphene resistance to deformations.The theory proves that graphene has two independent in-plane mechanical parameters and two independent out-of-plane mechanical parameters.Subsequently, Huang et al. obtained the deformation energy density of hexagonal boron nitride (h-BN) using the DREIDING force field, and also proved that the monolayer h-BN has four independent mechanical parameters [26].By combining the classical fracture theory and the interaction potential of carbon atoms, the researchers in [27] theoretically explained the brittle fracture of graphene.The above studies demonstrate that the macroscopic continuum mechanics theory needs modification to describe the mechanical behaviors of nanomaterials.
The existing MoS 2 molecular dynamics (MD) calculations have shown that the bending stiffness obtained using classical plate theory with the thickness of the three layers of atoms (h = 3.2 nm) is not identical to the stiffness obtained using MD calculations [28].To solve this contradiction, Huang proposed a nonlinear plate theory with independent in-plane and out-of-plane mechanical parameters to model MoS 2 mechanical behaviors based on finite temperature [29].This theory has a deformation energy density similar to classical plate theory, but it has four independent mechanical parameters.This new theory abandons the equivalent thickness of MoS 2 and directly takes in-plane and out-of-plane stiffnesses as independent mechanical parameters.Consequently, the Yakobson paradox is effectively avoided.
Two-dimensional nanomaterials are typically sensitive to temperature because their out-of-plane stiffness is low.MoS 2 expands with increasing temperature [30,31].Recent MD calculations have shown that temperature changes have little influence on the elastic parameters of MoS 2 ; however, temperature can cause significant thermal expansion [32].For single-layer MoS 2 with immovable boundaries, thermal expansion may induce thermal stress, which can lead to thermal buckling.In this study, the nonlinear static bending and vibrations of single-layer MoS 2 with four hinged edges were investigated based on a modified plate model proposed by Huang [29].This study focuses on the influence of temperature on the nonlinear mechanical behavior of monolayer MoS 2 .

Materials and Methods
A modified Föppl-von Karman plate model with independent in-plane and out-ofplane stiffnesses was established by Huang to model the mechanical properties of singlelayer MoS 2 [29].Because Huang's theory was published in Chinese, we briefly review this new theory for reader understanding.
Single-layer MoS 2 is considered a 2D plate in Huang's theory, as shown in Figure 1, and its deformation energy density is as follows [29]: ( ) ( ) Here, H and K are the mean and Gaussian curvatures, respectively, of the de- formed MoS2 middle surface.
( ) Using the von Karman nonlinear strain, the components of the strain tensor can be expressed as follows: Here, u , v , and w are the displacements of the middle surface in the x , y , and z directions, respectively.Equation ( 3) is consistent with the classical plate in addition to the four stiffness parameters.Therefore, we can define the Airy function of the in-plane 2D stress as Therefore, the in-plane strain is expressed as follows: Here, H and K are the mean and Gaussian curvatures, respectively, of the deformed MoS 2 middle surface.Q = det ε 0 ij and J = tr ε ij are the two invariants of the 2D strain tensor ε 0 ij , i, j = x, y on the middle surface.k B and k G are independent bending stiffness and torsional stiffness (Gaussian stiffness), respectively, whereas k b and k g are in-plane stiffness parameters.These four independent stiffness parameters are obtained through atomic simulations and experiments.
Using the von Karman nonlinear strain, the components of the strain tensor can be expressed as follows: Here, u, v, and w are the displacements of the middle surface in the x, y, and z directions, respectively.Equation ( 3) is consistent with the classical plate in addition to the four stiffness parameters.Therefore, we can define the Airy function of the in-plane 2D stress as N xx = ∂ 2 F/∂y 2 and N yy = ∂ 2 F/∂x 2 .From Equation (2), we obtain Therefore, the in-plane strain is expressed as follows: where According to Equation ( 5), the in-plane strain energy density can be rewritten as follows: To ensure the continuity and single-value of the displacement field, the strain field must satisfy the following completeness condition [33] This equation can be rewritten as ∆ 2 F = −χK.The Lagrange multiplier l(x, y) must be introduced into the potential energy function because the stress function is introduced.Then, Equation (2) can be rewritten as follows: By performing complex but direct computations on Equation (7) and identifying the Lagrange multiplier, Equation ( 7) can be transformed into Considering the influence of temperature on MoS 2 , Huang applied a boundary axial external force and thermal stress to the structure [29]; therefore, the load work is as follows: where q(x, y, t) is the load in the z direction, and N 0 xx and N 0 yy are the pre-applied axial tensile stresses in the x and y directions at the boundaries, respectively.N T xx and N T yy are the thermal stresses in the x and y directions at the boundaries [29,34] The Lagrange function can be constructed as L = U − W, which is subjected to variational calculations, namely, letting δL = 0 with the independent variables w and F. Therefore, we obtain the force balance equation and compatibility condition as follows: (10) Equation ( 10) is a mathematical model of MoS 2 derived by Huang [29], in which the out-of-plane and in-plane stiffness parameters are independent.To study the dynamics problem, we add the inertial force term m ∂ 2 w/∂t 2 to Equation (10) using the D'Alembert principle [35]; thus, the move equations can be rewritten as follows: For simplification, we assume that q(x, y, t) in Equation ( 11) is a harmonic load; therefore, q(x, y, t) = f (x, y) cos(Ωt).( 13) Here, we define the following dimensionless variables as follows: where a and b represent the side lengths of the monolayer MoS 2 , as shown in Figure 1.Thus, Equation (11) can be simplified into a dimensionless form as follows:

Analysis of Static Bending
Because Equations ( 15) and ( 16) are nonlinear partial differential equations, they are difficult to solve accurately.Therefore, the Galerkin method [36] is used to transform Equations ( 15) and ( 16) into ordinary differential equations in time.Equations ( 15) and ( 16) resemble the classical plate (but the mechanical parameters of the MoS 2 are independent).Under small deformations, symmetric loads may only induce symmetric deformations, even when nonlinear terms emerge.We then analyze the static and dynamic bending deformations under symmetric loads through first-and third-order symmetric modes.We expand the transverse displacement w and stress function F as follows: w = u 1 t sin π x sin π y + u 3 t sin 3π x sin π y, F = ξ 11 t sin π x sin π y + ξ 31 t sin 3π x sin π y. ( Substituting the F in Equation (17) into Equation ( 16) and multiplying sin π x sin π y and sin 3π x sin π y on the two sides of Equation ( 16) (the Galerkin model), we have The parameters in Equation ( 18) are as follows: Similarly, we substitute the w in Equation (17) into Equation (15); subsequently, we multiply sin π x sin π y and sin 3π x sin π y on the two sides of Equation (15), and considering Equation (18), the vibration equations for the first-order and third-order modes can be obtained as follows: ..
The parameters in Equation (20) are listed in Appendix A. By omitting the inertial terms in Equation (20), the static deformations of the midpoint of the monolayer MoS 2 are obtained as follows: We take the MoS 2 mechanical parameters in Table 1 as an example.Using the data in the table and Equation (21), the static bending amplitudes under external loads are obtained with four hinged edges, as shown in Figure 2a,b.
Table 1.Mechanical parameters of the single-layered MoS 2 [28,37,38].The two figures show that the thermal stress and axial compressive stress decrease the MoS2 stiffness.If both the first-and third-order modes are considered, the deflections of the midpoint would be slightly smaller than those considering only the first-order mode, as shown in Figure 2b.This may indicate that the first-order mode can precisely represent static deformations under symmetric loads.

Nonlinear Primary Resonance without Internal Resonance
In this section, our study of the nonlinear vibrations of single-layer MoS2 using Equation (20) are presented.These equations only contain cubic nonlinear terms.If the geometric dimensions of MoS2 and the axial force have the given values, the equations may exhibit a 1:3 internal resonance.This study mainly includes the following three parts: exciting only the first-or third-order primary resonance and the 1:3 internal resonance with the load's frequency near the low-order natural frequency.
To simplify the analysis, we assume that the damping force is 2 j j c u  with damping coefficient ˆj c .If there is no internal resonance in Equation ( 20), the vibrations of the unexcited modes rapidly decay because of damping.Therefore, the steady-state vibrations only contain the excited mode [39].Thus, the coupling terms in Equation ( 20) can be neglected.The forced vibration equation of single-layer MoS2 in the first-order mode ( where . We employ the multiple scale perturbation method to analyze Equation (22).The method is a classical perturbation method that is used to solve weak nonlinear differential equations.We supposed that the influences of damping, the nonlinear terms, and the loads emerge in a unified perturbation equation, so we set 22) is rewritten as follows: ( ) The two figures show that the thermal stress and axial compressive stress decrease the MoS 2 stiffness.If both the first-and third-order modes are considered, the deflections of the midpoint would be slightly smaller than those considering only the first-order mode, as shown in Figure 2b.This may indicate that the first-order mode can precisely represent static deformations under symmetric loads.

Nonlinear Primary Resonance without Internal Resonance
In this section, our study of the nonlinear vibrations of single-layer MoS 2 using Equation ( 20) are presented.These equations only contain cubic nonlinear terms.If the geometric dimensions of MoS 2 and the axial force have the given values, the equations may exhibit a 1:3 internal resonance.This study mainly includes the following three parts: exciting only the first-or third-order primary resonance and the 1:3 internal resonance with the load's frequency near the low-order natural frequency.
To simplify the analysis, we assume that the damping force is 2 ĉj .u j with damping coefficient ĉj .If there is no internal resonance in Equation ( 20), the vibrations of the unexcited modes rapidly decay because of damping.Therefore, the steady-state vibrations only contain the excited mode [39].Thus, the coupling terms in Equation ( 20) can be neglected.The forced vibration equation of single-layer MoS 2 in the first-order mode (Ω ≈ ω 1 ) or third-order mode (Ω ≈ ω 3 ) is as follows: .. where We employ the multiple scale perturbation method to analyze Equation (22).The method is a classical perturbation method that is used to solve weak nonlinear differential equations.We supposed that the influences of damping, the nonlinear terms, and the loads emerge in a unified perturbation equation, so we set ĉj = ε 2 c j and The small parameter ε = 0.1 is used in this research.Consequently, Equation ( 22) is rewritten as follows: ..
We assume the solution of Equation ( 23) as where T 0 = t, T 1 = εt, and T 2 = ε 2 t.By substituting Equation (24) into Equation ( 23) and equating the coefficients of ε and ε 3 on both sides, we obtain where In accordance with the ordinary differential equation theory, the solution of Equation ( 25) is as follows: where cc denotes the complex conjugate of the preceding term, and A j are the arbitrary functions of T 2 .Substituting Equation ( 27) into ( 26), we obtain where A is the complex conjugate of A. Letting Ω = ω j + ε 2 σ j and applying the elimination condition for the secular terms in Equation ( 28), we obtain We introduce the polar forms A j = (1/2)λ j exp iθ j with λ j and θ j in the real functions of T 2 .Substituting the A j into Equation (28) and separating the real and imaginary parts of Equation ( 29), we obtain where The steady-state motion will occur if D 1 λ j = D 1 γ j = 0.This corresponds to the solution of the following equations:

Primary Resonance of Low Frequency without Internal Resonance
The vibration amplitude of the low-frequency primary resonance can be obtained from Equation (31) using j = 1: One can obtain λ 1 from Equation (32).Then, we substitute it into Equation ( 24), so the first-order approximate solution is obtained.
The mechanical parameters of MoS 2 are listed in Table 1.The geometric dimensions are a = 5 nm and b = 10 nm, and the axial tensile stresses are N 0 xx = 0.3 nN/nm and N 0 yy = 0.1 nN/nm.These parameters prevent internal resonance.Because there has been no thorough research on damping, we use 2ε 2 c j = 0.05 for simplification.Therefore, we have c j = 2.5, ĉj = 0.025 for ε = 0.1.By substituting these parameters into the expressions in Appendix A, we obtain 3ω 1 ≪ ω 3 .
When these parameters are substituted into Equation ( 32), the amplitude-frequency response curves at three different temperatures were obtained with f 11 = 10, as shown in Figure 3a.This figure shows that temperature had an insignificant effect on the MoS 2 amplitude-frequency response curve.However, the combination of the load frequency and temperature had a significant effect on the load-amplitude curve, as shown in Figure 3b. Figure 3a,b also shows that the amplitude-frequency and load-amplitude response curves have two bifurcation points that lead to jumps in the vibration amplitude.The dotted lines in Figures 3 and 4 indicate unstable solutions.The stability of steady-state solutions can be determined through the eigenvalues of the Jacobian matrix of Equation (30); the details can be found in [39].
Materials 2024, 17, x FOR PEER REVIEW 9 of 21 ( ) The mechanical parameters of MoS2 are listed in Table 1.The geometric dimensions are  32), the amplitude-frequency response curves at three different temperatures were obtained with 11 10 f = , as shown in Figure 3a.This figure shows that temperature had an insignificant effect on the MoS2 amplitude-frequency response curve.However, the combination of the load frequency and temperature had a significant effect on the load-amplitude curve, as shown in Figure 3b. Figure 3a,b also shows that the amplitude-frequency and load-amplitude response curves have two bifurcation points that lead to jumps in the vibration amplitude.The dotted lines in Figures 3 and 4 indicate unstable solutions.The stability of steady-state solutions can be determined through the eigenvalues of the Jacobian matrix of Equation (30); the details can be found in [39].( ) The mechanical parameters of MoS2 are listed in Table 1.The geometric dimensions are  32), the amplitude-frequency response curves at three different temperatures were obtained with 11 10 f = , as shown in Figure 3a.This figure shows that temperature had an insignificant effect on the MoS2 amplitude-frequency response curve.However, the combination of the load frequency and temperature had a significant effect on the load-amplitude curve, as shown in Figure 3b. Figure 3a,b also shows that the amplitude-frequency and load-amplitude response curves have two bifurcation points that lead to jumps in the vibration amplitude.The dotted lines in Figures 3 and 4 indicate unstable solutions.The stability of steady-state solutions can be determined through the eigenvalues of the Jacobian matrix of Equation (30); the details can be found in [39].To validate the effectiveness of the approximate analytical solution, we simulate Equation ( 22) with f 11 , N 0 xx , N 0 yy = (10, 0.3, 0.1) using the Runge-Kutta method at T = 0 and T = 40, as shown in Figure 4a,b.Comparing the approximate analytical and numerical solutions, we find that the approximate analytical solution has good accuracy.The unstable solutions are indicated by dotted lines in Figure 3a,b and Figure 4a.

Primary Resonance of High Frequency without Internal Resonance
For the primary high-frequency resonance, the vibration amplitude can be obtained from Equation (31) using j = 3.We square the two equations and add them such that Substituting the λ 3 and γ 3 determined by Equation (34) into Equation (24), we obtain the third-order approximate solution as follows: The frequency-response or load-amplitude curves can be obtained using Equation (34) at different temperatures when the third-order mode is excited, as shown in Figure 5a,b.The precision of the approximate analytical solution is examined using the Runge-Kutta method as shown in Figure 6a,b.The dotted lines in Figures 5 and 6 indicate the unstable solutions.The stability of the steady-state solutions can be determined using the eigenvalues of the Jacobian matrix of Equation (30); the details can be found in [39].
From Figure 5a,b and Figure 6a, the following three main results can be drawn.First, the vibration amplitudes of the third-order mode are significantly smaller than those of the first-order mode under the same load.Second, the temperature has little effect on the amplitude of the third-order mode.Finally, the bifurcation points of the f 11 of the lowfrequency primary resonance are much smaller than those of the high-frequency resonance.A small bifurcation point of the first-order mode indicates that the low-frequency vibration is more prone to a large vibration amplitude.We use the Runge-Kutta method to calculate Equation (22) with f 31 = 50, as shown in Figure 6a      Here, we have used two thicknesses [28], h = 0.445 nm and h = 0.65 nm, to show the differences in nonlinear vibrations between the classical Föppl-von Karman plate model and the modified Föppl-von Karman plate model in this paper.So, the partial mechanical parameters with different effective thicknesses are shown in Table 2.The outstanding differences in the frequency-response curves between the two models can be found from Figure 7a,b.

Primary Resonance and 1:3 Internal Resonance at Low Frequency
To research the 1:3 internal resonance, we rewrite Equation ( 20) as follows:

Primary Resonance and 1:3 Internal Resonance at Low Frequency
To research the 1:3 internal resonance, we rewrite Equation ( 20) as follows: Here, we add two damping terms to the first-and third-order equations.The solution of Equation ( 36) are as follows: By substituting Equation (37) into Equation ( 36) and equating the coefficients of ε and ε 3 , we obtain ε : According to the theory of ordinary differential equations, we assume that the solution of Equation ( 38) are as follows: Here, cc denotes the complex conjugate of the preceding terms, and A 1 A 2 represent the functions of T 2 .Substituting Equation (40) into (39), we obtain where A denotes the complex conjugate of A and NST is the non-secular term.These equations may exhibit an internal resonance of ω 3 ≈ 3ω 1 .Introducing the detuning parameters σ 1 and σ 2 , we obtain (43) Therefore, the solvability conditions of Equations ( 41) and (42) are as follows: If the polar coordinate form A m = (1/2)a m exp(iβ m ), m = 1, 2 is introduced and substituted into Equation (45), and the real and imaginary parts of the equation are separated, we obtain Here, α m and β m are real functions of T 2 , and The steady-state motion may occur when D 2 a m = D 2 γ m = 0.The steady-state solution can be obtained using the following nonlinear equations: An algebraic equation revealing the relationship between the vibration amplitude and other parameters can be derived by squaring the second and fourth formulas of Equation ( 47) and summing them as follows: Letting a 2 1 = A 1 , Equation ( 48) can be simplified as a cubic equation for A 1 , where (50) Equation ( 49) can be rewritten as follows: with A 1 = x − e/3d, p = 3dg − e 2 /3d 2 , and q = 27d 2 h −9deg + 2e 3 /27d 3 .According to the Cardano formula, the solutions of Equation ( 51) are as follows: (52) There is one real root and two complex roots in Equation (51) if ∆ ≜ (q/2) 2 + (p/3) 3 > 0. When ∆ = 0 and p, q ̸ = 0, there is one double root and one single root; the equation has three distinct real roots if ∆ < 0.
A 1 is real because it is the vibration amplitude.Thus, we disregard the complex roots.Equation (48) provides the relationship between the low-frequency vibration amplitude a 1 and the high-frequency vibration amplitude a 2 .Equation (52) implies that a given a 2 corresponds to one or three values of a 1 .To simplify the research, we only consider that an a 2 corresponds to an a 1 .The more intricate cases will be studied in another paper.Hence, we had Equation ( 53) can be transformed into From the second and the fourth formulas of Equation ( 46), we obtain Substituting sin γ 1 and cos γ 1 into the first and third formulas of Equation (47), and then squaring and summing the two equations, we obtain The vibration amplitude of the first-and third-order models can be obtained using the following procedure: First, a 2 1 is obtained using Equation (54); subsequently, a 2 1 is substituted into Equation (54) to find a 2 .The angles γ 1 and γ 2 can be obtained by substituting a 1 and a 2 into the second and fourth formulas of Equation (47).

Stability Analysis of Steady-State Solutions
The stability of the solutions can be determined by investigating the nature of the singular points in Equation (46).To accomplish this, we set Subsequently, by substituting them into Equation (46) and considering that α 0 j , γ 0 j , j = 1, 2 meet Equation (47), we obtain The Jacobian matrix of Equation ( 57) is as follows: The elements in this matrix are shown in Appendix B. From Equation (47), we obtain ± π, The stability of steady-state solutions can be determined using the eigenvalues of the Jacobian matrix J.The steady-state solution is unstable if the eigenvalues of the corresponding Jacobian matrix contain positive real components.The unstable solutions are indicated by the dotted lines in Figures 8-11.
To study the 1:3 internal resonance, we use the mechanical and geometric coefficients listed in Tables 1 and 3; the axial stresses are N 0 xx = 10 nN/nm, N 0 yy = 49.8 nN/nm.By substituting these data into Equations ( 54) and ( 56), the relationship between the load frequency and amplitude can be displayed with a 1:3 internal resonance, as shown in Figures 8-11.third-order modes with three external forces at T 0 = .The two figures show that the vi- bration amplitudes of the low-order mode are significantly larger than those of the highorder mode under the same loads.Because the nonlinear terms increase the stiffness of MoS2, the resonance peaks shift toward higher frequencies.Furthermore, the vibration amplitude may increase with frequency, which indicates significant changes in motion.Figure 9a,b shows the effects of temperature on the amplitude-frequency response curves with 11 100 f = .The two figures reveal that the temperature has a more significant influence on the higher-order mode than on the lower-order mode.A slight temperature difference may induce an abrupt increase in the vibration amplitude.third-order modes with three external forces at T 0 = .The two figures show that the vi- bration amplitudes of the low-order mode are significantly larger than those of the highorder mode under the same loads.Because the nonlinear terms increase the stiffness of MoS2, the resonance peaks shift toward higher frequencies.Furthermore, the vibration amplitude may increase with frequency, which indicates significant changes in motion.Figure 9a,b shows the effects of temperature on the amplitude-frequency response curves with 11 100 f = .The two figures reveal that the temperature has a more significant influence on the higher-order mode than on the lower-order mode.A slight temperature difference may induce an abrupt increase in the vibration amplitude.To validate the reliability of the approximate analytical solutions, we perform numerical calculations for Equation (20) using the Runge-Kutta method.A comparison between the analytical solution and numerical calculations is shown in Figure 10a (the firstorder mode) and Figure 10b (the third-order mode).The results indicate that the approximate analytical solution is reliable.
To show the effect of temperature on the vibration amplitude, we draw the loadamplitude response curves with 2 7 σ = under different temperatures, as shown in Figure 11a,b.They imply that temperature has a more significant impact on the higher-order mode than on the low-order mode.The time-history curves display this feature, as shown in Figure 12a,b at 11 100 f = .
- To validate the reliability of the approximate analytical solutions, we perform numerical calculations for Equation (20) using the Runge-Kutta method.A comparison between the analytical solution and numerical calculations is shown in Figure 10a (the firstorder mode) and Figure 10b (the third-order mode).The results indicate that the approximate analytical solution is reliable.
To show the effect of temperature on the vibration amplitude, we draw the loadamplitude response curves with 2 7 σ = under different temperatures, as shown in Figure 11a,b.They imply that temperature has a more significant impact on the higher-order mode than on the low-order mode.The time-history curves display this feature, as shown in Figure 12a,b at 11 100 f = .
- Figure 8a,b illustrates the amplitude-frequency response curves of the first-and thirdorder modes with three external forces at T = 0.The two figures show that the vibration amplitudes of the low-order mode are significantly larger than those of the high-order mode under the same loads.Because the nonlinear terms increase the stiffness of MoS 2 , the resonance peaks shift toward higher frequencies.Furthermore, the vibration amplitude may increase with frequency, which indicates significant changes in motion.
Figure 9a,b shows the effects of temperature on the amplitude-frequency response curves with f 11 = 100.The two figures reveal that the temperature has a more significant influence on the higher-order mode than on the lower-order mode.A slight temperature difference may induce an abrupt increase in the vibration amplitude.
To validate the reliability of the approximate analytical solutions, we perform numerical calculations for Equation (20) using the Runge-Kutta method.A comparison between the analytical solution and numerical calculations is shown in Figure 10a (the first-order mode) and Figure 10b (the third-order mode).The results indicate that the approximate analytical solution is reliable.
To show the effect of temperature on the vibration amplitude, we draw the loadamplitude response curves with σ 2 = 7 under different temperatures, as shown in Figure 11a,b.They imply that temperature has a more significant impact on the higher-order mode than on the low-order mode.The time-history curves display this feature, as shown in Figure 12a

Conclusions
In this study, we employ a modified plate model in which four independent elastic parameters and thermal stresses are considered to investigate the nonlinear static bending and vibrations of monolayer MoS2.First, we use the Galerkin method to truncate the partial differential equation with the first and third modes.Subsequently, nonlinear static bending and forced vibrations are explored using ordinary differential equations obtained using the Galerkin method.The main conclusions are as follows: (1) The first-order mode can accurately represent the static deformation of MoS2 under symmetric loads.(2) Temperature has a slight effect on the single-mode vibrations of the MoS2.However, the combination of the load frequency and temperature have a more significant effect on the vibrations.When the temperature has a slight change, the bifurcation points of vibration amplitude will change significantly with the identical load's amplitude and frequency.(3) The bifurcation points of the load at the low-frequency primary resonance are significantly smaller than those at the high frequency for single-mode vibrations.(4) The vibration amplitudes of the first-order mode are significantly larger than those of the higher-order modes under the same loads when a 1:3 internal resonance appear in the MoS2.(5) For the 1:3 internal resonance, the temperature has a more significant influence on the higher-order mode than on the lower-order mode, and a slight temperature difference may induce an unexpected jump in the vibration amplitude.Under the same load, the maximum value of the amplitude-frequency curve will increase significantly with the temperature's increase.
The above findings may give some important inspirations when a single-layer MoS2 is used in nano-resonators and mass sensors.

Conclusions
In this study, we employ a modified plate model in which four independent elastic parameters and thermal stresses are considered to investigate the nonlinear static bending and vibrations of monolayer MoS 2 .First, we use the Galerkin method to truncate the partial differential equation with the first and third modes.Subsequently, nonlinear static bending and forced vibrations are explored using ordinary differential equations obtained using the Galerkin method.The main conclusions are as follows: (1) The first-order mode can accurately represent the static deformation of MoS 2 under symmetric loads.(2) Temperature has a slight effect on the single-mode vibrations of the MoS 2 .However, the combination of the load frequency and temperature have a more significant effect on the vibrations.When the temperature has a slight change, the bifurcation points of vibration amplitude will change significantly with the identical load's amplitude and frequency.(3) The bifurcation points of the load at the low-frequency primary resonance are significantly smaller than those at the high frequency for single-mode vibrations.(4) The vibration amplitudes of the first-order mode are significantly larger than those of the higher-order modes under the same loads when a 1:3 internal resonance appear in the MoS 2 .(5) For the 1:3 internal resonance, the temperature has a more significant influence on the higher-order mode than on the lower-order mode, and a slight temperature difference may induce an unexpected jump in the vibration amplitude.Under the same load, the maximum value of the amplitude-frequency curve will increase significantly with the temperature's increase.

Appendix B
The Jacobian matrix coefficients of Equation (58) are as follows:

εFigure 1 .
Figure 1.Calculation diagram of a single-layer MoS2 under load: (a) plate model with the coordinate; (b) applied edge loads; (c) side view of the MoS2 lattice structure.

Figure 1 .
Figure 1.Calculation diagram of a single-layer MoS 2 under load: (a) plate model with the coordinate; (b) applied edge loads; (c) side view of the MoS 2 lattice structure.
. The thermal stresses in the uniform temperature field are N T xx = k b ε T xx = k b αT and N T yy = k g ε T yy = k g αT, where ε T xx and ε T yy represent the thermal strain, and α is the coefficient of thermal expansion (CTE).

.
These parameters prevent internal resonance.Because there has been no thorough research on damping, we use ε = .By substituting these parameters into the expres- sions in Appendix A, we obtain When these parameters are substituted into Equation (

Figure 3 .Figure 4 .
Figure 3. (a) Frequency-response curves of low-frequency primary resonance with the three temperatures for 11 10 f = ; (b) load-response curves of vibration amplitudes with two temperatures.

1 .Figure 3 .
Figure 3. (a) Frequency-response curves of low-frequency primary resonance with the three temperatures for f 11 = 10; (b) load-response curves of vibration amplitudes with two temperatures.

.
These parameters prevent internal resonance.Because there has been no thorough research on damping, we use When these parameters are substituted into Equation (

Figure 3 .Figure 4 .
Figure 3. (a) Frequency-response curves of low-frequency primary resonance with the three temperatures for 11 10 f = ; (b) load-response curves of vibration amplitudes with two temperatures.

Figure 5 .
Figure 5. (a) Frequency-response curves of the high-frequency primary resonance with three temperatures for 31 50 f = ; (b) load-response curves of vibration amplitudes with two temperatures.

3 Figure 5 .
Figure 5. (a) Frequency-response curves of the high-frequency primary resonance with three temperatures for f 31 = 50; (b) load-response curves of vibration amplitudes with two temperatures.

Figure 5 .Figure 6 .
Figure 5. (a) Frequency-response curves of the high-frequency primary resonance with three temperatures for 31 50 f = ; (b) load-response curves of vibration amplitudes with two temperatures.

Figure 7 .
Figure 7. (a) Frequency-response curves of the low-frequency primary resonance with the classical plate model and the modified plate model" for 10 f = ; (b) frequency-response curves of the high- frequency primary resonance with the classical plate model and the modified plate model for 50 f = .

Figure 5 .Figure 6 .
Figure 5. (a) Frequency-response curves of the high-frequency primary resonance with three temperatures for 31 50 f = ; (b) load-response curves of vibration amplitudes with two temperatures.

Figure 7 .
Figure 7. (a) Frequency-response curves of the low-frequency primary resonance with the classical plate model and the modified plate model" for 10 f = ; (b) frequency-response curves of the high- frequency primary resonance with the classical plate model and the modified plate model for 50 f = .

λ 1 σ 1 Figure 7 .
Figure 7. (a) Frequency-response curves of the low-frequency primary resonance with the classical plate model and the modified plate model" for f = 10; (b) frequency-response curves of the highfrequency primary resonance with the classical plate model and the modified plate model for f = 50.

Figure
Figure 8a,b illustrates the amplitude-frequency response curves of the first-and

Figure
Figure 8a,b illustrates the amplitude-frequency response curves of the first-and

Figure 10 .Figure 11 .
Figure 10.(a) Comparison between approximate analytical and numerical solutions for vibration amplitudes of the first-order model with two loads for T 0 = ; (b) comparison between approxi- mate analytical and numerical solutions for vibration amplitudes of the first-order model with two temperatures for 11 100 f = .

Figure 10 .Figure 10 .Figure 11 .
Figure 10.(a) Comparison between approximate analytical and numerical solutions for vibration amplitudes of the first-order model with two loads for T = 0; (b) comparison between approximate analytical and numerical solutions for vibration amplitudes of the first-order model with two temperatures for f 11 = 100.

Figure 11 .Table 3 .
Figure 11.(a) Amplitude-response curves of the first-order mode for three temperatures at σ 2 = 7; (b) amplitude-response curves corresponding to the third-order mode for three temperatures at σ 2 = 7.

Figure 12 .
Figure 12.(a) Time-history response curves of the first-order mode with two temperatures at 11 100 f = ; (b) time-history response curves of the third-order mode with two temperatures at 11 100 f = .