Investigation of critical load of structures using modified energy method in nonlinear-geometry solid mechanics problems

2.1 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 .

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:

2.2 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, U = { u , v , w } T , intra-element displacements ( U e ) 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 .

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.

2.3 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.

2.3.1 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:

k ( i − 1 ) t + Δ t Δ D ( i ) = F ext t + Δ t − F int ( i − 1 ) t + Δ t ,

(12) D ( i ) t + Δ t = D ( i − 1 ) t + Δ t + Δ D ( i ) .

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 :

k ( 0 ) t + Δ t = k t ,

F int ( 0 ) t + Δ t = F int t ,

(13) D ( 0 ) t + Δ t = D t .

Using the following equation, define a residual or unbalanced force, F R , at each iteration:

(14) F R ( i ) t + Δ t = F ext t + Δ t − F int ( i − 1 ) t + Δ t .

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.

Figure 1

Equilibrium equations solved through original (top) and modified Newton–Raphson methods (bottom).

2.3.2 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.

Figure 2

Response of structures with a damage point.

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:

(15) F int t + Δ t = λ t + Δ t F ext .

Defining the internal load at the i th iteration of an iterative process as follows:

(16) F int ( i ) t + Δ t = F int ( i − 1 ) t + Δ t + k ( i − 1 ) t + Δ t Δ D ( i ) ,

Also, the loading factor is as follows:

(17) λ ( i ) t + Δ t = λ ( i − 1 ) t + Δ t + Δ λ ( i ) .

For these methodologies, the primary load equilibrium equation takes the iterative form as shown in Eq. (18):

(18) k ( i − 1 ) t + Δ t Δ D ( i ) = ( λ ( i − 1 ) t + Δ t + Δ λ ( i ) ) F ext − F int ( i − 1 ) t + Δ t .

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):

(19) f ( Δ D ( i ) , Δ λ ( i ) ) = 0 .

Now, assuming that displacement, δ D ( i ) , and load factor, δ λ ( i ) , vary per load step to the i th iteration as follows:

(20) δ D ( i ) = D ( i ) t + Δ t − D t , δ λ ( i ) = λ ( i ) t + Δ t − λ t .

Then, as illustrated in Figure 3, the following constraint [33] can be defined by using a constant spherical arc length, Δ ℓ :

(21) ( δ λ ( i ) ) 2 + ( δ D ( i ) ) T ( δ D ( i ) ) β = Δ ℓ 2 .

Figure 3

Geometrical description of 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:

(22) Δ ℓ Min ≤ Δ ℓ ≤ Δ ℓ Max .

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:

(23) Δ D ¯ ( i ) = ( k ( i − 1 ) t + Δ t ) − 1 ( λ ( i − 1 ) t + Δ t F ext − F int ( i − 1 ) t + Δ t )  , Δ D ¯ ¯ ( i ) = ( k ( i − 1 ) t + Δ t ) − 1 ( F ext ) .

The displacement increment at the i th iteration can be expressed as follows using these two equations:

(24) Δ D ( i ) = Δ D ¯ ( i ) + Δ λ ( i ) Δ D ¯ ¯ ( i ) .

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. (25).

(25) Δ λ ( i ) = − ( F ext ) T Δ D ¯ ( i ) ( F ext ) T Δ D ¯ ¯ ( i ) .

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.

4.1 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 E A = 3 × 10 5 N 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.

Figure 5

Powell’s tower [36].

Figure 6

Load–displacement curve at Powell’s tower’s Node 1 [36].

Figure 7

Comparing the CL of Powell’s tower using the Powell’s and the energy methods (MEM).

4.2 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², while members 22–35 had a cross-sectional area of A _22–35 = 51.61 cm². 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].

Figure 8

Thompson–Hunt truss structure [37].

Figure 9

Load–displacement curve at Thompson–Hunt truss structure’s Node 10.

Figure 10

Comparison of CL at Node 10 of the Thompson–Hunt truss structure by MEM, Torkamani and Shieh [39], and Kondoh and Atluri [38].

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].

Figure 11

Load–displacement curve at Thompson–Hunt truss structures’ Node 19.

Figure 12

Comparison of CL at Node 19 of the Thompson–Hunt truss structure by MEM, Torkamani and Shieh [39], and Kondoh and Atluri [38].

4.3 Four-member truss structure

In this example, the effect of vertical load on a four-member 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 schematically depicts this truss structure in conjunction with its applied loads. E A = 9 × 10 4 N 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.

Figure 13

Four-member truss structure [40].

Figure 14

Vertical load–displacement curve at Node 2 of a four-member truss structure via the TLA and MEM [40].

Figure 15

TLA and MEM comparison of the CL of a four-member truss structure.

4.4 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 E A = 3 × 10 5 N 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.

Figure 16

Eight-member truss structure [41].

Figure 17

Horizontal load–displacement curve at Node 1 along the eight-member truss structure [41].

Figure 18

Comparison of CL at Node 1 of the eight-member truss structure using the Narayanan and MEM methods.

Figure 19

Horizontal load–displacement curve at Node 7 along the eight-member truss structure [40].

Figure 20

Comparison of CL at Node 7 of the eight-member truss structure using the Narayanan and MEM methods.

4.5 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.

Figure 21

Arched truss structure [38].

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.

Figure 22

Vertical load–displacement curve at Node 10 along the arched truss structure.

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].

Greek

β

contribution factor of velocity to displacement

δ

increment

Δ

delta of the quadratic equation of velocity/change in the corresponding parameter

Δ ℓ

arc length

Δ ℓ Min

upper bound to arc length

Δ ℓ Max

lower bound to arc length

Δ t

timestep

ε e

intra-element strain

λ

loading factor

σ e

intra-element stress

Latin

A

the cross-section area of members

A , B , C

coefficients of quadratic equation

B

strain–displacement relation matrix

b

the horizontal length of members

D

displacement vector of the system

d e

nodal displacement vectors of element

E

Lagrangian strains

E

elasticity modulus

E int

the internal energy of structure

E ext

load-induced energy

e

Eulerian strains

F int

internal forces vector

F ext

external forces vector

F R

imbalanced forces vector

f S

resistant force

h

truss height

I

moment of inertia of members

k t

tangent stiffness

N

shape functions vector

p

point load

t

pseudo-time

U

displacement vector

V e

element volume

v ( t )

velocity at the time t

x

current position of system

X

the initial position of the system