Investigation of critical load of structures using modi ﬁ ed energy method in nonlinear - geometry solid mechanics problems

: Geometrically nonlinear analysis is required for resolving issues such as loading causes failure and struc - ture buckling analysis. Although numerical methods are recommended for estimating the exact solution, they lack the necessary convergence in the presence of bifurcation points, making it challenging to ﬁ nd the equilibrium path using these methods. Thus, the modi ﬁ ed energy method is employed instead of the numerical method, frequently used to solve quasi - static problems with nonlinear nature and bifurcation points. The ultimate goal of this study is to determine the critical load of structures through the modi ﬁ ed energy method rather than other methods in which the relationship between force, displacement, and constraint is used to solve the problem. This study ﬁ rst describes the energy method for this type of problem and then details its computational steps progressively. This method yields numerical results when applied to numerical examples such as truss and frame structures and coded in MATLAB software. These ﬁ ndings are com - pared to the analytical results. The energy method is more precise than the alternative methods and superior to the Newton – Raphson method at crossing the load – displace - ment curve ’ s bifurcation points.


Introduction
Numerous physical phenomena occur in nature that is nonlinear. Under certain conditions, conventional linear models can approximate nonlinear models. A common assumption made in the linearization of solid mechanics problems is that displacements or strains are smallan assumption that produces results that are adequate for only a few practical applications. Indeed, the primary advantage of this assumption is that it maintains the solution of the linearized problem to be nearly identical to that of the nonlinear analysis problem, even though solving the linearized problem is fundamentally easier. However, when the structure is subjected to large deformations, linearization no longer yields adequate results, and the problem must be solved in terms of the nonlinearity encountered in reality. For instance, when a structure is subjected to an earthquake event that causes it to deviate from its linear region (i.e., significant structural deformations), it is no longer reasonable to assume constant structural stiffness and nonlinear interaction between applied load and structural deformations. Finding a suitable solution is only impossible in this case when a nonlinear analysis is used. In general, because finding an exact solution is impractical for most nonlinear problems, numerical techniques can be viewed as viable alternatives to analytic methods for such problems.
Various researchers have examined numerous nonlinear problems involving related structures. Generally, the leading research on nonlinear analysis of structures has been dedicated to the geometrically nonlinear analysis of buckling [1]. For instance, in several early works, an incremental technique was used to combine a geometric stiffness matrix and a method for updating the coordinate system [2,3]. The nonlinear effects of materials were investigated for the first time in refs. [4,5], with the concept of tangent stiffness matrix introduced in refs. [5,6]. Notably, until that time, all solution procedures relied heavily on the forward-Euler method, resulting in extremely high calculation error levels. The Newton-Raphson method was invented several years later to address this latter issue. Later on, a modified Newton-Raphson method was proposed to improve the efficiency of the original Newton-Raphson method by eliminating the need for frequent updating of the stiffness matrix [4,7]. The direct energy search was presented in the same vein of research [8,9]. Finally, Zienkiewicz proposed using subincrements to improve the solution accuracy [10].
Massicotte and Fafard (1991) investigated the geometric analysis application of the arc-length method. As a result, they developed two numerical methods using continuous and discrete formulations of the arc-length method, with the disadvantages and advantages of each method analyzed further. According to the findings of this study, some of Crisfield's method's drawbacks could be addressed to accelerate convergence to the solution [11]. Coda and Greso (2004) discussed several studies involving a straightforward formulation for large deformations in two-dimensional (2D) frames. Their work focused on the position as the primary unknown parameter, and their formulation was entirely based on a spatial description of the body while remaining independent of displacement concepts. Strains were calculated in their study by estimating the difference between the original and deformed forms of a structure in a one-dimensional (1D) domain, using the strain energy at all points as a starting point.
The case studies revealed that this method's primary advantage is its simplicity and ability to analyze nonlinear frames in a 2D domain [12]. Rabczuk et al. (2006) introduced a method for analyzing shear strips¹ with an adhesive surface. They based their method on the planar strain model, emphasizing the material's nonlinear effects on adhesive surfaces (e.g., concrete), and applied it in 2D and three-dimensional (3D) domains for materials with rate-dependent and rate-independent deformation. Finally, the results indicated that the proposed method could perform nonlinear analysis in shear strip modeling (for interpreting the material type and analyzing the structural geometry) [13]. Lorenzana et al. (2011) proposed a simple yet effective methodology for the nonlinear analysis of 2D frames.
Nonlinear effects were assumed to be concentrated at the beam's ends, with the remainder behaving linearly in this formulation. They then developed an elastoplastic model to explain the phenomenon of failure. Finally, a planar frame was subjected to nonlinear analysis by defining a liquefaction function (load distribution) in terms of axial load, shear force, and flexural momentum. They sought to determine the ultimate load on 2D frames by considering a failure model concentrated at the beam's two ends [14].
Mansouri and Safari proposed a method similar to the Newton-Raphson method but considered geometrical nonlinearities. This method, commonly referred to as the two-point iteration method, begins by generating an initial guess and then iterating the solution. To this end, a mathematical function is selected to approximate the tangent matrix. In other words, this methodology combines the Newton-Raphson and two-point procedures in order to reduce the processing time and iterations. Following the elimination of trial-and-error steps, the nonlinear analysis of planar frames demonstrated a 40% reduction in the processing time compared to the Newton-Raphson method [15].
Numerous researchers have published on the subject of this study more recently. Mamouri et al., for example, investigated nonlinear geometric instabilities in 2D frames. This study assumes that nonlinear instabilities in elastic materials are associated with large deformations prior to buckling. Total Lagrange equations were used to solve the problem. This methodology is based on a concept termed dynamic relaxation, which allows the discovery of equilibrium points beyond the limiting points, a significant advantage over conventional static methods [16]. Radnic et al. (2016) compared several numerical models with 1D and 2D elements for nonlinear static analysis of concrete frames. Their research focused on the nonlinear analysis of planar frames, and the equations considered were independently written for steel and frame, with the analytical results and experimental data compared further [17]. Moita et al. (2016) proposed a nonlinear model for compact plates to analyze the sandwich structures' nonlinear behavior. Under the assumption of continuous displacement along boundaries, layers of defined thickness were used as elements in the finite-element method (FEM), with the layers having distinct functions. Classical plates modeled stiff elastic layers, while Reddy's third-order shear deformation theory simulated the core. The equilibrium path was then determined using the Newton-Raphson method with frequent increments. Additionally, the arc-length method was used to determine the correct path when the load-deformation curve bifurcates.
Overall, various reports have indicated that the proposed model could more effectively reduce the processing load than a 3D finite-element analysis by utilizing flat triangular elements [18]. The authors of ref. [19] conducted  nonlinear analyses on cabled structures and proposed a hybrid formulation. This research introduced a series of chain elements based on a consistent yet imprecise the combination of finite deformations and finite-element modeling and included a continuous axial force element and a distributed axial force element. The curvilinear coordinate system was used to investigate the problem at hand in this contribution. Additionally, the principle of virtual work was used to derive the relatively weak relationship between strain and displacement. Finally, the FEM was examined through the presentation of three examples. In comparison to the previous equations that assumed a constant axial force, the proposed equations with variable axial force have the potential to yield accurate results.
Feng et al. published an analysis of Timoshenko's bending-shear beam with fiber elements in 2017. This study used a softening model with a plasticity component to simulate failure. With the superposition of shear and bending in mind, a beam element based on Timoshenko's beam was selected and then investigated using the displacement function. The material behavior was modeled using the fiber section method, in which two distinct types of sections were formed using steel and concrete fibers. Finally, it was determined that the effect of superposition on the shear behavior of beams could be adequately investigated using the elements of this beam through analysis and comparison of experimental results on simple beams [20]. Vafaei et al. (2019) established the compact plasticity model as a reliable predictor of the nonlinear behavior of reinforced concrete frames subjected to vertically increasing loads [21].
The authors of ref. [22] proposed a numerical procedure for developing a methodology for geometrically nonlinear problems encountered in solid mechanics. The emphasis was primarily on the transition from load-based equilibrium equations to step-by-step energy equilibrium procedures. In most cases, bifurcation points are unavoidable in nonlinear problems. As a result, this methodology placed a premium on the method's accuracy and convergence, given its numerical nature. The proposed methodology's outputs were compared to exact solutions (if any) and other methods' outputs to ensure validity. Notably, the investigation was conducted under the assumption of quasi-static loading, which ignored the effects of damping and inertia. Moreover, the forces were susceptible to a change in direction due to structural deformationsa change that can impose nonlinear effects. Nonetheless, the latter nonlinear effect was ignored, and the structure was assumed not to affect boundary conditions. Razaghi et al. (2019) investigated the numerical efficiency of the energy method for predicting the prebuckling response of framed structures with bifurcation points [23]. Mohemmi et al. (2020) proposed an equivalent method based on the results of pullout tests for simulating rod slippage in reinforced concrete frames [24]. For practical analyses, the results of FEM analysis in Abaqus suggested a steel strength reduction factor of 0.6. Li et al. (2019) developed a model for analyzing the plastic collapse of reinforced concrete frames with shear walls [25]. This study modeled the fracture process using the cohesive crack concept and the traction-separation law. Lezgy-Nazargah et al. developed a refined global-local FEM for bending analysis of laminated composite rectangular beams in another study [26]. The presented element was free of shear locking and satisfied the continuity conditions for displacements, and transverse shear stresses at interfaces, as well as the nonhomogeneous shear tractions on the beam's upper and lower surfaces. By comparing the results of this study to those of the previous research and the 3D theory of elasticity, it was determined that this formulation is accurate for thin and relatively thick beams.
Lezgy-Nazargah (2017) developed a time-consuming and effective 1D FEM to assess the progressive failure of laminated composite beams based on the refined highorder global-local beam theory [27]. To predict the firstply failure load, they used an iterative method. Following that, a stiffness degradation factor was used to modify the material properties of the failed elements. Finally, the applied load was increased incrementally for each load increment until the laminate reached its ultimate strength. Another study developed a materially nonlinear model for reinforced concrete beams.
In contrast to the commonly used Euler-Bernoulli or Timoshenko kinematic assumptions for 1D beam elements, a layered global-local theory was used to account for shear deformation effects normal flexibility. The nonlinear behavior of cracked concrete and reinforcing bars was also simulated using fixed smeared crack and elastoplastic models. Overall, a comparison of the results revealed that the proposed method produced accurate results while requiring minimal computational resources [28].
Einafshar et al. (2021) presented a FEM to analyze the beam's flexural, axial buckling, buckling, and nonlinear geometric properties. Nonclassical effects such as cross-section and conventional flexibility are incorporated into the formulation via a novel structural concept termed equivalent layer cross-sectional modeling. A layer of the composite section replaces the initial section of the thin-walled beam with equivalent stiffness in this framework. The displacement fields of the beam are represented using the layer shear deformation theory. The Newton-Raphson linearization method is used to solve nonlinear equations. The proposed method does not require the use of the shear correction factor and has a low degree of freedom (DOF). Comparing the results demonstrates that the proposed formula for stability and nonlinear geometric analysis is precise and efficient [29].
Lezgy-Nazargah et al. (2019) presented a FEM for analyzing insulated concrete composite foam sandwich panels. Nonlinear structures can be predicted using 1D elements. The numerical results demonstrate that the 1D elements used in this study can accurately calculate the flexural capacity and stress distribution for all loading stages [30].

Techniques for solving quasistatic nonlinear problems
This section begins with a brief overview of nonlinear problems encountered in solid mechanics and then examines several common practical techniques for resolving such problems.

Different types of nonlinear problems in solid mechanics
Generally, the nonlinearity exhibited by solid materials can be classified under two categories: material nonlinearity or geometric nonlinearity. Stress levels are not proportional to strains in the presence of material nonlinearity. Geometric nonlinearity occurs when the deformation of a structural system increases to the point where the deformed system is fundamentally different from the original. The equilibrium equation can no longer be expressed regarding the structure's initial geometry. A structural system can exhibit four distinct types of nonlinearity. When writing the cinematic relationship between system displacement and corresponding strains, the nonlinearity is geometric if a nonlinear term can be introduced. If the relationship between stress and strain is not linear, the nonlinearity is material in nature. Practitioners rarely consider two additional nonlinearities concerning boundaries and applied forces. Strains cannot be calculated as a linear function of displacements in geometrically nonlinear problems; otherwise, unrealistically large deformations for rigid body rotation may be calculated, or nonunique strains may be measured. However, it is important to note that large displacements do not always imply large strains across the structure. For example, the analysis of flexible thin-walled structures could be considered. In general, when performing structural analysis in the presence of geometric nonlinearity, one of two descriptions, total Lagrange or updated Lagrange, can be used. The total Lagrange is suitable for problems involving large displacements and small strains, whereas the updated Lagrange is appropriate when the system encounters large (plastic) strains [31].
Significant differences between the initial geometry and the deformed structure geometry occur due to applied loads in a gas of geometric nonlinearity across a system. This is in contrast to linear analyses, which calculate all stresses, strains, and displacements based on the system's initial geometry. Because the structure is continuously deformed while remaining in equilibrium, equilibrium equations must be based on the system's deformed geometry. In general, when deformations are large (geometrically nonlinear analysis), the difference between the initial and current configurations (i.e., deformed geometry) is too large to ignore. As a result, understanding how to describe these large deformations, stresses, and strains in a material is critical in such cases. Consider a solid body in this case; large displacement results in a change in geometry from the initial state (nondeformed geometry) to the current configuration (deformed geometry) [31]. All formulas and relationships in this study are denoted by bold vector or matrix parameters.
The displacement vector, u, is defined in Eq. (1) as a function of the system's current position, x, and the initial system geometry, X.

= − u x X.
(1) Given the definition of strain as the derivative of displacement concerning an axis, both the initial and formed geometries can be used in the derivation, resulting in two distinct approaches to strain definition. If the initial geometry is used as a reference point for the derivation, the Lagrangian strains, E, are defined as in Eq. (2): Alternatively, Eulerian strains, e, can be expressed as follows, using the initial geometry as a reference for derivation:

Equilibrium equation in a nonlinear state
The virtual work principle, which is widely used in the formulation of nonlinear problems, states that when a structure is in equilibrium, the virtual works performed by internal and external forces are equal for each consistent virtual displacement to the system. This section reviews the procedure for determining the equilibrium equation for a nonlinear FEM. In a finite-element problem with a 3D displacement field, can be defined as shape functions, N , as shown in Eq. (4).
where d e is the element's nodal displacement vector. It can provide us with the intra-element displacements, ε e according to matrix Eq. (5), when using the definition of matrix B.

= ε
Bd . e e (5) By utilizing the variational principle and specifying a virtual field for the element's displacements in the form of δd e , Eq. (5) can be rewritten as Now, by equating the internal work (caused by stresses, σ e ) to the integration over the element volume, V e , and the external work (caused by nodal loads, F e ), the virtual work principle implies the following: The equilibrium equation can be written as follows after some matrix simplification: After obtaining the equilibrium equation at the element level, assembling was used to obtain the system's general equilibrium equations, as follows: where F int and F ext denote the system's internal and external forces, respectively.

Solving nonlinear equilibrium equations
Numerous numerical methods for solving the equilibrium equations of a nonlinear structural system have been proposed. The Newton-Raphson iterative method combined with an incremental-iterative solution strategy is a well-known traditional approach to solving such problems. However, some more advanced techniques have been proposed to improve the efficiency of solution techniques for various problems. Examples of such algorithms include linear searches, path-following, quasi-Newton, and branch switching techniques at bifurcation points.
As mentioned in the previous section 2.2, the equilibrium equations in solid mechanics relating internal and external forces applied to a solid body, i.e., Eq. (9), can be rewritten to describe the situation at time + t t Δ as follows: Given that the internal force is a nonlinear function of the system displacements, D, variations in the internal force, F Δ int , between times t and + t t Δ can be expressed using the concept of tangent stiffness coefficient, k t , as defined in Eq. (11): Previous experience indicates that if you proceed directly from the D Δ calculation to the next step, the solution will contain significant errors, or the approach will become unstable. As a result, it is necessary to proceed with the calculations in one step until an allowable tolerance is reached to ensure accurate solutions. The following section discusses several conventional methods in this field.

Newton-Raphson method
The Newton-Raphson method is a well-known numerical technique for the analysis of nonlinear systems and can be expressed as follows: The superscript i denotes the iteration number in the preceding equations, with the following conditions considered at the start of the iteration process, i.e., = i 1: Using the following equation, define a residual or unbalanced force, F R , at each iteration: A stopping condition for the iteration can be expressed as a lower bound on this out-of-equilibrium load vector to ensure that the calculations are completed in one step.
A consideration is that each iteration of the tangent stiffness matrix incurs substantial computational cost. By recognizing this limitation, a modified Newton-Raphson algorithm was proposed in which the tangent stiffness matrix was calculated only once at the start of each time step and then remained unchanged throughout the iterations at that time step. Figure 1 illustrates these two methodologies geometrically. Other techniques, commonly referred to as quasi-Newton methods, have been proposed for the same purpose. They represent transient approaches between the original and modified Newton-Raphson methods.

Load-displacement constraint methods
Nonlinear analyses are frequently used to determine the critical damage load applied to a structural system. In this case, as illustrated in Figure 2, large increments for the load can be considered in the initial portion of the load-displacement curve, where the structure exhibits linear behavior. As one approaches the damage point (the point at which the load-displacement relationship becomes nonlinear), the load increments must decrease, particularly when passing through the damage point (corresponding to zero stiffness), where mathematical singularity occurs. Additionally, the postbuckling response of the structure must be determined using special methods that can accommodate both load reduction and displacement increase concurrently.
Riks [32] pioneered the concept of load-displacement constraint methods. This methodology generally introduces a loading factor λ into the equilibrium equations to increase/decrease the applied force's intensity. Indeed, the inclusion of this factor in the equations accomplishes three major goals: (i) accelerating convergence at each time step, (ii) facilitating passage through the damage point, and (iii) allowing for the discovery of a postdamage solution. Equilibrium equations are rewritten as follows in this series of methods: Defining the internal load at the ith iteration of an iterative process as follows: Also, the loading factor is as follows:  For these methodologies, the primary load equilibrium equation takes the iterative form as shown in Eq. (18): Eq. (18) contains two significant unknowns: displacement increment, ( ) D Δ i , and load increment, ( ) Δλ i , which must be related to one another via a constraint, as illustrated in Eq. (19): Now, assuming that displacement, ( ) D δ i , and load factor, ( ) δλ i , vary per load step to the ith iteration as follows: Then, as illustrated in Figure 3, the following constraint [33] can be defined by using a constant spherical arc length, ℓ Δ : In the preceding equation, β is a constant used to create a dimensionless expression. In practice, an upper and a lower bounds ( ℓ Δ Max and ℓ Δ Min , respectively) on the arc length are imposed: When the structure exhibits near-linear behavior, values close to ℓ Δ Max are selected, whereas values near ℓ Δ Min must be adopted when the structural response enters a nonlinear zone.
Equilibrium equations are occasionally rewritten as follows: The displacement increment at the ith iteration can be expressed as follows using these two equations: Substituting these equations into the constraint shown in Eq. (21) results in a quadratic equation for ( ) λ Δ i . This equation must be appropriately initialized before the system can be analyzed. Additionally, if the constant increment in external work is used, the constraint mentioned above will take the form of Eq.
It is necessary to emphasize that the methodologies discussed in this section are merely initiators. This implies that a different solution algorithm should be used to start the solution process. Furthermore, additional techniques are required for the proper selection of ℓ Δ during the analysis process. Moreover, whenever the algorithm detects motion toward divergence, it must cease iteration before restarting the solution process with new parameters.

Research method
The methodology used in this study is based on the modified energy method (MEM)a technique first introduced in refs. [34,35] for the time history analysis of nonlinear structures. MEM was coined in general [34] for determining the dynamic response of structures with one DOF, including linear and simple nonlinear systems such as a spring with third-degree stiffness (Duffing), Coulomb friction, and pendulum motion at large rotations (nonlinear sinusoidal term). The adopted methodology begins by discretizing the energy equilibrium equations in time and then using force equilibrium equations to determine the system's true velocity.
Later, this methodology was applied to shear frame systems with linear behavior in ref. [35], where the elimination of discontinuous velocities was first introduced in the solution of coupled equations obtained using energy methods on the structures mentioned above. This concept (the application of an energy-based approach) is extended in the current research to analyze geometrically nonlinear problems encountered in solid mechanics (which include significantly more complicated nonlinear terms than those addressed in the latter two references). We generally compare MEM's performance to that of other methods (e.g., Newton-Raphson and Riks-Wampner) when dealing with problems involving bifurcation points along their load-displacement curve. As previously stated, a realistic analysis of a structural problem must include dynamic analysis of the structural behavior over time. However, in many cases, when the structure is loaded incrementally, a quasi-static process is assumed to read analytic complexities (as in this study). In a quasi-static analysis, the structure is loaded incrementally. In such cases, the variable t denotes a pseudo-time variable indicating the loading intensity at a specific time step rather than the conventional time variable used in nonlinear dynamic analyses.
Instead of using Eq. (9) to represent the system's force equilibrium, this study employed an energy-based formulation as the governing equation for a nonlinear structural problem. In general, the energies present in a static system are composed of the energy induced by the structural strength, F int , and the energy/work induced by the loading on the structure, F ext , which are denoted herein as E int and E ext , respectively. According to energy conservation and Eq. (26), these two energy components must be equal at all times.
The preceding equation can be rewritten as follows to obtain an incremental equation: where t j and + t j 1 represent two successive points in time. In general, energy changes induced by internal work, E Δ int , for a nonlinear system can be described in the form of the following integration: The above equation can be discretized as follows using the trapezoidal law, where j denotes the time step number: Eq. (29) takes the following form when the internal force at the current step, + F j int 1 , is expressed in terms of the internal force at the previous step, F j int , and the tangent stiffness matrix, k t : In contrast, the work performed by the external work, E Δ ext , can be described as Eq. (31): Using the definition of the velocity, v j , the above equation can be rewritten as follows: Which takes the following form upon discretization: In contrast, displacement changes, D Δ , can be expressed using the Eulerian equation, as demonstrated in Eq. (35), where β is the coefficient of the velocity contribution to determining the system displacements, i.e., determines the type of integration. as a function of + v j 1 could be obtained.
Coefficients of this equation are defined in the following steps: Given that the presence of imaginary velocities in a structural system implies that they have no physical meaning, the delta of the quadratic Eq.
, must always be positive. Given that this value is generally dependent on the system's general characteristics, applied load, and discretization parameters, a constraint that requires the numerical results to be stable may be added. This is a significant advantage that this formulation has over more traditional load-based methods.
As shown in Eq. ), thereby eliminating the possibility of analyzing the structural system with imaginary velocities. However, for any other value (i.e., < ≤ β 0 1), the delta sign should be monitored throughout the analysis, and the length of the time step t Δ reduced if the sign becomes negative. Indeed, Eq. (37) demonstrates unequivocally that the delta of a quadratic equation is always nonnegative as the time step approaches zero (i.e., → t Δ 0). Another question to address is how, given the quadratic nature of the instantaneous velocity equation, which is supposed to return two velocity values at each instant, the true system velocity can be selected from the two. According to ref. [30], when velocity is considered continuous, the true velocity can be chosen as the one closest to the velocity at the previous time step. To solve the quadratic equation of velocity and determine the true velocity of the system at each instant using this methodology, we assume continuous changes in velocity over time, thereby ruling out the physical consistency of discontinuous velocity changes. As a result, the velocity closest to the velocity in the preceding time step is chosen as the true velocity of the structure from the two velocities obtained from these quadratic equations. As a result, this technique is referred to in the literature to eliminate discontinuous velocities. By denoting the calculated velocities in the current time step as v 1 and v 2 , we obtain the following: Generally, the following steps should be followed when analyzing a nonlinear system for computer execution: (Figure 4) The MEM has several advantages over other methods for calculating the critical load (CL) of structures: -Indirect integration methods necessitate the selection of regulatory parameters. The values for these parameters that are recommended are based on a simplified hypothesis. When dealing with nonlinear thematic problems, it is debatable whether to select and calibrate the constants in this type of method.
-Along with their highly complicated mathematical relationships, methods based on the principle of energy equilibrium have several limitations, including the requirement for stability, the absence of nonlinear sentences, and, most importantly, their limitation to structures with a single DOF. Additionally, they are applicable only when there are almost no nonlinear constraints.
-None of the methods of nonlinear analysis' optimal time value selection provide the analyst with any information about obtaining an accurate and consistent result. Moreover, the recommended values for t derived from various procedures are based solely on linear analyses and do not consider the effect of system load type and analysis time on the response stability method.
Appropriate t Δ selection in nonlinear analysis, providing information about the structure's stability, and finally calculating the structure's CL are also unique characteristics of this numerical method not found in other practice methods. Finally, by employing an energy-based formulation and discontinuous velocity elimination, it is intended to derive a structure (with general behavior) from a nonlinear analysis without calculating the length of the nonlinear analysis.

Numerical examples
This section evaluated the proposed algorithm from the previous section's efficiency on several numerical Solving the quadraƟc equaƟon of velocity using Eq.36 Defining the actual speed of the system using Eq.38 CalculaƟng displacements in the current step ( j+1 D= j D+∆D)

UpdaƟng structural sƟffness ( t K)
Have the answers examples by coding in MATLAB. In the first instance, the unstable behavior of Powell's cantilever tower was investigated in the presence of large deformations and geometric nonlinearity. The Thompson-Hunt truss structure was considered in the following example, particularly in terms of the effect of imperfections on the truss's buckling. In the third example, a four-member truss structure with complicated geometric nonlinearity was considered. We examine an eight-member structure in a fourth example due to its highly nonlinear behavior. Finally, an arched truss was analyzed to determine the effect of geometrical imperfections. Of note, the proposed energy method's results were compared to those obtained using other numerical methods reported by other researchers and any exact analytic solutions for each example.

Powell's tower
The cantilevered tower proposed by Powell et al. was analyzed in this example [36]. This structure is forced into an unstable state by a combination of horizontal and vertical loads. The geometric configuration of this tower and the applied loads are depicted in Figure 5. The stiffness of the various members of this truss structure was assumed to be = × EA 3 10 N 5 in this section. After analyzing this tower, the load-displacement curve for Node 1 is shown in Figure 6, indicating that the energy method is highly accurate, as demonstrated by the excellent agreement between the proposed methodology's results and those reported in ref. [36] using high-order strains. The CLs obtained using the two methods are shown in Figure 7.

Thompson-Hunt truss structure
Thompson and Hunt initially analyzed this structure to determine the effect of possible imperfections on the buckling of 3D grid trusses and trussed columns [23]. Kondoh and Atluri investigated the relationship between local buckling and the structure's overall stability [38].
The Thompson and Hunt truss structure with 35 members and 19 nodes is depicted in Figure 8. This structure's members all have circular cross-sections and have an elasticity modulus of E = 68.964 GPa. Members 1-21 had a cross-sectional area of A 1-21 = 54.84 cm 2 , while members 22-35 had a cross-sectional area of A 22-35 = 51.61 cm 2 . This test determines the overall stability of the structure when a horizontal compressive load is applied at Node 19. Of note, this structure is not symmetrical horizontally, with its neutral axis located somewhere above the structure's horizontal center. Vertical displacement at Node 10 and horizontal displacement at Node 19 were calculated after analyzing the structure. The vertical load-displacement plot for Node 10 and the horizontal load-displacement curve for Node 19 are shown in Figures 9 and 10, respectively. These figures illustrate the energy method's results and other researchers' results [38,39].  A closer examination of these figures demonstrates that the proposed methodology produces extremely precise results when analyzing geometric nonlinearity. The CL of the structure is compared (i.e., validated) in Figures 11  and 12 using the proposed energy method and procedures introduced by Torkamani and Shieh [39] and Kondoh and Atluri [38].

Four-member truss structure
In this example, the effect of vertical load on a fourmember truss structure is analyzed [40]. One rigid and one soft support were used in this truss structure. This structure exhibits complex geometric nonlinearity as a result of its support configuration. Figure 13 [39], and Kondoh and Atluri [38].
depicts this truss structure in conjunction with its applied loads. = × EA 9 10 N 4 was the member stiffness. Vertical displacement at Node 2 was calculated through structural analysis. The vertical load-displacement curve at Node 14 is illustrated in Figure 14 using the energy method and the total Lagrangian formulation (usually known as the total Lagrangian approach [TLA]). This figure demonstrates an acceptable agreement between the proposed analysis method's results and those obtained using the total Lagrangian method with high-order Green strains. Additionally, Figure 15 illustrates the CL of the structure using the energy method and the total Lagrangian formulation, demonstrating the proposed energy method's acceptable accuracy.

Eight-member truss structure
An eight-member truss structure is analyzed in this example. Figure 16 illustrates the geometric configurations of this truss structure. The stiffness of the member was set to = × EA 3 10 N 5 in this case. Due to this truss structure's highly nonlinear behavior, numerous researchers have used it to evaluate their proposed formulations. A different variant of the same structure was analyzed by some researchers. Following this research, a methodology for analyzing this truss structure is developed, and the results are compared to those reported in [40]. Horizontal load-displacement curves for Nodes 1 and 7 were created for this purpose (indicated on the structure). The horizontal load-displacement curves for Nodes 1 and 7 are shown in Figures 17 and 18, respectively. Examining these figures reveals the proposed energy method's high accuracy for analyzing geometric nonlinearity, as demonstrated by the results' acceptable agreement with those of ref. [40]. Figures 19 and 20 compare the CLs obtained using the energy and Narayanan methods.

Arched truss structure
The arched truss structure with 35 members investigated in this study is depicted in Figure 21. Rosen and Schmit originally developed it to study the effects of geometric imperfections on structural behavior [42,43]. Kondoh and Atluri investigated the behavior of this structure in terms of local and total imperfections [38]. All truss members in this structure have circular cross-sections and an elasticity modulus of E = 68.964 GPa. The coordinates of all 19 truss nodes and the cross-sectional area of the structural members are listed in Tables 1 and 2. In this example, the overall behavior of the structure under a vertical compressive load was investigated at Node 10 by neglecting local and total imperfections. Vertical displacement at Node 10 was determined and plotted in Figure 22 following a structural analysis. The energy method's results are presented alongside those of other researchers for validation purposes.
A careful examination of these plots demonstrates that the proposed procedure accurately analyzes geometric  nonlinearities. Additionally, Figure 23 compares the CL values obtained using the energy method to those obtained using other procedures, illustrating the proposed methodology's high accuracy.       Figure 23: Comparison of vertical CL at Node 10 of the Thompson-Hunt truss structure by MEM, Torkamani and Shieh [39], and Kondoh and Atluri [38].

Conclusion
This study aimed to determine the numerical efficiency of the MEM for solving common quasi-static problems in solid mechanics. The energy method was used to solve numerical examples with geometric nonlinearity, and the following conclusions were reached: • Bifurcation points could be addressed using this method by correctly identifying the path to equilibrium on loaddeformation diagrams without encountering divergence when the structure has no stiffness. This is in contrast to Newtonian methods, which are incapable of performing in this manner. • The advantage of this method over the load-displacement constraint method is that it solves the problem by selecting the values of t Δ and β. This is because the presented quadratic equation allows for the observation of previous deltas. In contrast, the constraintbased method employs a trial-and-error approach to the same process. For example, in the arc-length method, the analysis must select ℓ Δ Min , ℓ Δ Max , β, and λ before the bifurcation points can be kept under control.
• Finally, a small number of sample structures that other researchers had previously analyzed were solved using the energy method, and the results were compared. This comparative study demonstrated the energy method's acceptable accuracy, establishing it as a promising method for structural analysis.  velocity at the time t x current position of system X the initial position of the system Acknowledgment: All the authors are grateful to the reviewers for their fruitful suggestions.
Funding information: The authors state that no funding was involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.