Generalized finite element formulation for efficient first-order plastic hinge analysis

Abstract

This article presents generalized finite element formulation for plastic hinge modeling based on lumped plasticity in the classical Euler–Bernoulli beam. In this approach, the plastic hinges are modeled using a special enrichment function, which can describe the weak discontinuity of the solution at the location of the plastic hinge. Furthermore, it is also possible to insert a plastic hinge at an arbitrary location of the element without modifying its connectivity or adding more elements. Instead, the formations of the plastic hinges are achieved by hierarchically adding more degrees of freedom to existing elements. Due to these features, the proposed methodology can efficiently perform the first-order plastic hinge analysis of large-frame structures. A generalized finite element solution technique based on the static condensation scheme is also proposed in order to reduce the computational cost of a series of linear elastic problems, which is in general the most time-consuming portion of the first-order plastic hinge analysis. The effectiveness and accuracy of the proposed method are verified by analyzing several representative numerical examples.

Introduction

It is well known that not a small number of structural collapse accidents have been occurring all over the world, and their main causes are natural and man-made disasters such as earthquakes, tsunamis, hurricanes, and terrorist attacks. A large number of civilian casualties and serious damage to properties are resulted in, as already have been observed in, for examples, Alfred P. Murrah federal building collapse in the United States in 1995, ¹ Taiwan earthquake in 1999, ² and Chile earthquake in 2010. ³ Therefore, it is very important to minimize the structural damages by these types of disasters and design structures that can effectively resist excessive loadings induced by them.

The use of numerical analysis methodologies such as the finite element method (FEM) is a very effective way to achieve these goals, and their importance has been increasing with the help of the recent rapid development of cutting-edge information technologies. It is possible to utilize massively parallel computing machines for obtaining accurate simulation results. However, it is more important to develop an efficient numerical methodology that can accurately simulate the collapse of building structures by consuming a relatively small amount of computational cost.

In general, the collapse of building structures is simulated by the successive formation of plastic hinges in the frame structure. Figure 1 illustrates an example of the progressive collapse of a building structure due to column loss, which frequently occurs in the structural collapse caused by explosion. ⁴ As another example, Figure 2 shows the progressive collapse of a frame structure subjected to earthquake excitation. ⁵ In this case, a large amount of internal bending moment is concentrated at the lower level of the structure, and it may lead to the sidesway failure mechanism as shown in the figure.

Figure 1.

Progressive collapse of a building structure due to column loss.

Figure 2.

Progressive collapse of a building structure due to earthquake excitation.

The plastic hinges can be modeled based on either lumped or distributed plasticity. The lumped plasticity model ideally assumes that plastic deformation is concentrated only at the location of plastic hinges and thus the plastic hinges can be described as an internal pin. The deformed shape of a frame member is continuous at the plastic hinge, but its derivative becomes discontinuous. The limit capacity of the plastic hinge is kept constant in the classical limit analysis, but it is possible to take into account the so-called softening plastic hinge concept with post-peak structural response as discussed in the work of Armero and Ehrlich,^6,7 which introduces the embedded discontinuity with local dissipative mechanism. This approach was subsequently applied to perform the collapse simulation of reinforced concrete ⁸ and steel frames.^9,10

In the distributed plasticity model, the member cross-section is composed of a number of fibers with axially nonlinear material properties, and the plastic deformation is spread over a certain length of the member near the plastic hinge, which is called a yield zone. ¹¹ As a result, there exists a strong gradient in the deformed shape of the member, but no discontinuity. Although the lumped plasticity model has limitation in accurately describing plastic deformation within a member, it is much simpler than the distributed plasticity model and can predict the overall collapse failure of a frame structure in the form of the first-order plastic hinge analysis. ¹²

This study presents a generalized finite element formulation for plastic hinge modeling based on lumped plasticity in the classical Euler–Bernoulli beam. In this approach, the plastic hinges are modeled using a special enrichment function, which can describe the weak discontinuity of the solution at the location of the plastic hinge. Furthermore, it is also possible to insert a plastic hinge at an arbitrary location of the element without modifying its connectivity or adding more elements. Instead, the formations of the plastic hinges are achieved by hierarchically adding more degrees of freedom (DOFs) to existing elements. Due to these features, the proposed methodology can efficiently perform the first-order plastic hinge analysis of large-frame structures. A generalized finite element solution technique based on the static condensation scheme is also proposed in order to reduce the computational cost of a series of linear elastic problems, which is in general the most time-consuming portion of the first-order plastic hinge analysis. This technique is similar to the one proposed by Gupta et al. ¹³ The effectiveness and accuracy of the proposed method are verified by analyzing several representative numerical examples.

The outline of this article is as follows. Section “Plastic hinge modeling using the GFEM” discusses the plastic hinge modeling based on lumped plasticity using the generalized finite element method (GFEM). A proper enrichment function is proposed to describe the weak discontinuity, which exists at the location of the plastic hinges. The effectiveness and accuracy of the proposed GFEM formulation are investigated by performing convergence analyses of a model problem using both the standard FEM and GFEM. In section “Efficient first-order plastic hinge analysis based on the GFEM,” the first-order plastic hinge analysis based on the GFEM is developed. This approach utilizes the static condensation scheme due to the hierarchical nature of generalized finite element shape functions, resulting in a significant amount of saving in the computational cost. Its effectiveness and accuracy are investigated by performing a first-order plastic hinge analysis on a model problem using both the standard FEM and GFEM. Section “Conclusion” provides the main conclusions of this study, and related future research topics are also suggested.

Plastic hinge modeling using the GFEM

In this section, the GFEM discretization for plastic hinge modeling is proposed, and its effectiveness is verified by performing numerical experiments and analyzing their results. As the numerical experiments, h- and p-extensions are performed on a model problem using both the standard FEM and proposed GFEM, and their convergence rates are compared with those of theoretical predictions.

Generalized finite element discretization for plastic hinge modeling based on lumped plasticity

This section discusses the plastic hinge modeling using GFEM shape functions and its characteristics. The classical Euler–Bernoulli beam is utilized for the first-order plastic hinge analysis performed in this study. The weak form of its governing equation is given by

$$\begin{array}{c}\mathrm{Find}{u}_h\in X(\Omega )\subset H^2(\Omega ),\Omega =\{x:0<x<H\}\hfill \\ \multicolumn{1}{c}{\mathrm{such}\mathrm{that},\forall \delta u_h\in X(\Omega )}\\ \underset{0}{\overset{H}{\int }}\mathrm{EI}\frac{{\mathrm{d}}^2{u}_h}{\mathrm{d}{x}^2}\frac{{\mathrm{d}}^2\delta u_h}{\mathrm{d}{x}^2}\mathrm{d}x=\underset{0}{\overset{H}{\int }}q\delta u_h\mathrm{d}x-{\left[\mathrm{EI}\frac{{\mathrm{d}}^3{u}_h}{\mathrm{d}{x}^3}\delta u_h\right]}_0^H\hfill \\ \multicolumn{1}{c}{+{\left[\mathrm{EI}\frac{{\mathrm{d}}^2{u}_h}{\mathrm{d}{x}^2}\frac{{\mathrm{d}}^2\delta u_h}{\mathrm{d}{x}^2}\right]}_0^H}\end{array}$$

(1)

where $X(\Omega )$ is a discretization of $H^2(\Omega )$ , a Hilbert space defined on $\Omega$ , built with generalized finite element shape functions. In addition, q is the uniformly distributed load, and H is the length of the beam. In the above equation, u_h is a generalized finite element approximation given by

$$\begin{array}{c}{u}_h\left(x\right)=\sum _{\alpha =1}^N{\phi }_{\alpha }\sum _{i\in \ell \left(\alpha \right)}(L_{\alpha i}{\underset{¯}{u}}_{\alpha i}^0+B_{\alpha }{L}_{\alpha i}{\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}})\hfill \\ =\underset{\mathrm{smooth}\mathrm{interpolation}}{\underbrace{\sum _{\alpha =1}^N\sum _{i\in \ell \left(\alpha \right)}{\phi }_{\alpha }{L}_{\alpha i}{\underset{¯}{u}}_{\alpha i}^0}}+\underset{\mathrm{non-smooth}\mathrm{interpolation}}{\underbrace{\sum _{\alpha =1}^N\sum _{i\in \ell \left(\alpha \right)}{\phi }_{\alpha }{B}_{\alpha }{L}_{\alpha i}{\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}}}}\hfill \end{array}$$

(2)

where ${\phi }_{\alpha }$ , $\alpha =1,\dots ,N$ are partition of unity functions, $L_{\alpha i},i\in \ell (\alpha )$ are polynomial enrichment functions defined at node α. In this equation, B_α is the plastic hinge enrichment defined at node α and required only at the nodes of elements with plastic hinges, and ${\underset{¯}{u}}_{\alpha i}^0$ and ${\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}}$ are the nodal DOFs corresponding to the polynomial and plastic hinge enrichments, respectively. The same shape functions are also used for the virtual displacement (δu_h) in the Galerkin method as follows

$$\delta u_h\left(x\right)=\sum _{\alpha =1}^N{\phi }_{\alpha }\sum _{i\in \ell \left(\alpha \right)}(L_{\alpha i}\delta {\underset{¯}{w}}_{\alpha i}^0+B_{\alpha }{L}_{\alpha i}\delta {\underset{¯}{w}}_{\alpha i}^E)$$

(3)

where $\delta {\underset{¯}{w}}_{\alpha i}^0$ and $\delta {\underset{¯}{w}}_{\alpha i}^E$ are the nodal DOFs of the virtual displacement approximation. Further details on the generalized finite element can be found in many references.^14–19

In the lumped plasticity approach, the displacement of the beam is continuous at the location of the plastic hinge, but its derivative is discontinuous, as illustrated in Figure 3. This behavior can be accurately described by the GFEM approximation expressed by equation (2). It consists of both the smooth and non-smooth interpolations. The smooth interpolation can represent the continuous deformation of the beam, while the non-smooth interpolation represents the discontinuity in the derivative of displacement (weak discontinuity) at the plastic hinge location. As can be seen from the second line of equation (2), the former is achieved by the polynomial enrichment functions and latter by the combination of the polynomial and plastic hinge enrichment functions.

Figure 3.

Plastic hinge description based on the lumped plasticity.

Since equation (1) contains the integral term of second derivative of displacement and virtual displacement requiring the C¹ feature of shape functions, we use the cubic Hermite shape functions associated with displacement as the partition of unity function ${\phi }_{\alpha }$ . They are the simplest partition of unity satisfying the C¹ feature as shown in Figure 4(a) and also enjoy the so-called Kronecker-delta property, that is

$${\phi }_{\alpha }\left(x_{\beta }\right)={\delta }_{\alpha \beta },(a,b=1,\dots ,N)$$

(4)

Figure 4.

GFEM shape functions for plastic hinge modeling: (a) partition of unity function and (b) enrichment for plastic hinge modeling.

In addition, the derivative of the partition of unity function defined at node α becomes zero at any node β, that is

$$\frac{\mathrm{d}{\phi }_{\alpha }\left(x_{\beta }\right)}{\mathrm{d}x}=0,(a,b=1,\dots ,N)$$

(5)

As the polynomial enrichment functions $(L_{\alpha i})$ , the following monomials are used

$$\begin{array}{c}{L}_{\alpha i}=\{1,\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right),{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^2\},(p=2)\hfill \\ \{1,\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right),{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^2,{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^3\},(p=3)\hfill \\ \{1,\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right),{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^2,{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^3,{\left(\frac{x-x_{\alpha }}{h_{\alpha }}\right)}^4\},(p=4)\hfill \end{array}$$

(6)

where $x_{\alpha }$ is the coordinate of node α, and $h_{\alpha }$ is a scaling factor. With the use of this form of the polynomial enrichment functions, the DOFs corresponding to polynomial enrichments with i = 1 and 2 ( ${\underset{¯}{w}}_{\alpha 1}^0$ and ${\underset{¯}{w}}_{\alpha 2}^0$ ) can have the physical meaning of displacement and rotation at node α, respectively, and thus the imposition of essential boundary conditions can become as straightforward as in the standard FEM. The normalization using $h_{\alpha }$ in equation (6) is helpful to minimize round-off errors, as already discussed by Duarte et al. ¹⁵

Figure 4(b) plots the plastic hinge enrichment function (B_α) included in equation (2), which is defined by

$$B_{\alpha }=\sum _{\alpha =1}^2{\phi }_{\alpha }\left|d_{\alpha }\right|-\sum _{\alpha =1}^2\left|{\phi }_{\alpha }{d}_{\alpha }\right|$$

(7)

where d_α is the distance between the plastic hinge location and node α. Only the nodes of the element having the plastic hinge, which are indicated by β and γ in the figure, are enriched with this special enrichment function. This function was originally proposed and used to describe the weak discontinuity at the bi-material interface by Moës et al., ²⁰ and it can also used for the description of the plastic hinge. This is further discussed in section “Implementation issues.” Notice that it satisfies C¹ requirement at the plastic hinge location and C^∞ at all the other locations. Furthermore, it can satisfy the following conditions, which can be effectively utilized for straightforward imposition of displacement boundary conditions

$$B_{\alpha }\left(x_{\beta }\right)=0,(a,b=1,\dots ,N)$$

(8)

$$\frac{\mathrm{d}{B}_{\alpha }\left(x_{\beta }\right)}{\mathrm{d}x}=0,(a,b=1,\dots ,N)$$

(9)

It can be also noted that, in general, any GFEM shape functions can be hierarchically added to an existing GFEM discretization as discussed by Duarte et al. ¹⁵ For example, in order to create a quartic GFEM approximation, it is possible to add only GFEM shape function terms associated with quartic monomials to the GFEM discretization constructed using linear to cubic GFEM approximations. Due to this feature of the method, the first-order plastic hinge analysis can be efficiently performed with the aid of the static condensation scheme introduced in section “Generalized finite element solution technique based on the static condensation.”

By plugging equations (2) and (3) into equation (1) based on the concept of the Galerkin method, we can obtain

$$\left(\delta {\mathbf{u}}^{\mathrm{T}}\right)\left\{\underset{-1}{\overset{1}{\int }}\mathrm{EI}{\mathbf{B}}^{\mathrm{T}}\mathbf{B}\frac{H}{2}\mathrm{d}\xi \right\}\mathbf{u}=\left(\delta {\mathbf{u}}^{\mathrm{T}}\right)\underset{-1}{\overset{1}{\int }}q{\mathbf{N}}^{\mathrm{T}}\frac{H}{2}\mathrm{d}\xi$$

(10)

where u and $\mathbf{B}$ can be defined from the generalized finite element approximation given by equation (2) as follows

$$\begin{array}{c}\multicolumn{1}{c}{u_h=\underset{\mathrm{smooth}\mathrm{interpolation}}{\underbrace{\sum _{\alpha =1}^N\sum _{i\in \ell \left(\alpha \right)}{\phi }_{\alpha }{L}_{\alpha i}{\underset{¯}{u}}_{\alpha i}^0}}+\underset{\mathrm{non}-\mathrm{smooth}\mathrm{interpolation}}{\underbrace{\sum _{\alpha =1}^N\sum _{i\in \ell \left(\alpha \right)}{\phi }_{\alpha }{B}_{\alpha }{L}_{\alpha i}{\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}}}}}\\ \multicolumn{1}{c}{=\left[\begin{array}{cc}{\mathbf{N}}^0\hfill & {\mathbf{N}}^{\mathrm{ph}}\hfill \end{array}\right]\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]=\mathbf{Nu}}\end{array}$$

(11)

and

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\mathrm{d}}^2{u}_h}{\mathrm{d}{x}^2}=\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\mathbf{Nu}=\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\left[\begin{array}{cc}{\mathbf{N}}^0\hfill & {\mathbf{N}}^{\mathrm{ph}}\hfill \end{array}\right]\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]}\\ \multicolumn{1}{c}{=\left[\begin{array}{cc}{\mathbf{B}}^0\hfill & {\mathbf{B}}^{\mathrm{ph}}\hfill \end{array}\right]\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]=\mathbf{Bu}}\end{array}$$

(12)

Here, ${\underset{¯}{\tilde{\mathbf{u}}}}^0$ is the column matrix of DOFs associated with polynomial enrichments, ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ is that associated with plastic hinge enrichments, and ${\mathbf{N}}^0$ and ${\mathbf{N}}^{\mathrm{ph}}$ are the GFEM shape function matrices associated with ${\underset{¯}{\tilde{\mathbf{u}}}}^0$ and ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ , respectively. Similarly, ${\mathbf{B}}^0$ and ${\mathbf{B}}^{\mathrm{ph}}$ are the second derivatives of ${\mathbf{N}}^0$ and ${\mathbf{N}}^{\mathrm{ph}}$ , respectively. As in the case of the standard FEM, a linear mapping between the natural $(\xi )$ and global (x) coordinate systems is adopted to facilitate numerical integration of the weak formulation. Finally, equation (10) can be simplified into

$$\left[\begin{array}{cc}{\mathbf{K}}^0\hfill & {\mathbf{K}}^{0,\mathrm{ph}}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}^{\mathrm{ph},0}}& {\mathbf{K}}^{\mathrm{ph}}\hfill \end{array}\right]\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}^0\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}^{\mathrm{ph}}}\end{array}\right]$$

(13)

where ${\mathbf{K}}^0={\int }_{-1}^{1}EI{({\mathbf{B}}^0)}^{\mathrm{T}}{\mathbf{B}}^0(H/2)\mathrm{d}\xi$ , ${\mathbf{K}}^{0,ph}={\int }_{-1}^{1}EI{({\mathbf{B}}^0)}^{\mathrm{T}}{\mathbf{B}}^{ph}(H/2)\mathrm{d}\xi$ , ${\mathbf{K}}^{\mathrm{ph},0}={\int }_{-1}^1\mathrm{EI}{({\mathbf{B}}^{\mathrm{ph}})}^{\mathrm{T}}{\mathbf{B}}^0(H/2)\mathrm{d}\xi$ , and ${\mathbf{K}}^{ph}={\int }_{-1}^{1}EI{({\mathbf{B}}^{ph})}^{\mathrm{T}}{\mathbf{B}}^{ph}(H/2)\mathrm{d}\xi$ . Similarly, ${\mathbf{F}}^0={\int }_{-1}^{1}q{({\mathbf{N}}^0)}^{\mathrm{T}}(H/2)\mathrm{d}\xi$ , and ${\mathbf{F}}^{\mathrm{ph}}={\int }_{-1}^{1}q{({\mathbf{N}}^{\mathrm{ph}})}^{\mathrm{T}}(H/2)\mathrm{d}\xi$ .

Implementation issues

This section discusses two issues related to the implementation of the GFEM approximation proposed in section “Generalized finite element discretization for plastic hinge modeling based on lumped plasticity,” which are the imposition of displacement boundary conditions and numerical integration of the proposed GFEM discretization.

Imposition of displacement boundary conditions

Another advantage of the scaled monomials expressed in equation (6) is that the nodal DOFs corresponding to polynomial enrichments with i = 1 and 2 ( ${\underset{¯}{u}}_{\alpha 1}^0$ and ${\underset{¯}{u}}_{\alpha 2}^0$ ) can have a physical meaning of the displacement and rotation at node α, respectively. For example, if the scaled monomials up to p = 3 are used as the enrichment functions at all nodes of the domain, due to equations (4) and (8), the GFEM approximation of displacement at node β can be computed by

$$\begin{array}{c}\multicolumn{1}{c}{u_h\left(x_{\beta }\right)=\sum _{\alpha =1}^N{\phi }_{\alpha }\left(x_{\beta }\right)\sum _{i\in \ell \left(\alpha \right)}(L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^0+B_{\alpha }\left(x_{\beta }\right)L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}})}\\ ={\phi }_{\beta }\left(x_{\beta }\right)\sum _{i\in \ell \left(\beta \right)}(L_{\beta i}\left(x_{\beta }\right){\underset{¯}{u}}_{\beta i}^0+B_{\beta }\left(x_{\beta }\right)L_{\beta i}\left(x_{\beta }\right){\underset{¯}{u}}_{\beta i}^{\mathrm{ph}})\hfill \\ =1·[1·{\underset{¯}{u}}_{\beta 1}^0+\left(\frac{x_{\beta }-x_{\beta }}{h_{\beta }}\right)·{\underset{¯}{u}}_{\beta 2}^0+{\left(\frac{x_{\beta }-x_{\beta }}{h_{\beta }}\right)}^2·{\underset{¯}{u}}_{\beta 2}^0\hfill \\ \multicolumn{1}{c}{+{\left(\frac{x_{\beta }-x_{\beta }}{h_{\beta }}\right)}^3·{\underset{¯}{u}}_3^0]}\\ ={\underset{¯}{u}}_{\beta 1}^0\hfill \end{array}$$

(14)

This indicates that ${\underset{¯}{u}}_{\alpha 1}^0$ has the physical meaning of the vertical displacement at node β. Similarly, due to equations (4), (5), (8), and (9), the approximated rotation of the beam can be computed by

$$\begin{array}{c}{\theta }_h\left(x_{\beta }\right)=\frac{\mathrm{d}{u}_h\left(x_{\beta }\right)}{\mathrm{d}x}\hfill \\ =\sum _{\alpha =1}^N\frac{\mathrm{d}{\phi }_{\alpha }\left(x_{\beta }\right)}{\mathrm{d}x}\sum _{i\in \ell \left(\alpha \right)}(L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^0+B_{\alpha }\left(x_{\beta }\right)L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}})\hfill \\ +\sum _{\alpha =1}^N{\phi }_{\alpha }\left(x_{\beta }\right)\sum _{i\in \ell \left(\alpha \right)}\frac{\mathrm{d}}{\mathrm{d}x}(L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^0+B_{\alpha }\left(x_{\beta }\right)L_{\alpha i}\left(x_{\beta }\right){\underset{¯}{u}}_{\alpha i}^{\mathrm{ph}})\hfill \\ ={\phi }_{\beta }\left(x_{\beta }\right)\sum _{i\in \ell \left(\beta \right)}\hfill \\ \multicolumn{1}{c}{(\frac{\mathrm{d}{L}_{\beta i}\left(x_{\beta }\right)}{\mathrm{d}x}{\underset{¯}{u}}_{\beta i}^0+\frac{\mathrm{d}{B}_{\beta }\left(x_{\beta }\right)}{\mathrm{d}x}{L}_{\beta i}\left(x_{\beta }\right){\underset{¯}{u}}_{\beta i}^{\mathrm{ph}}+B_{\beta }\left(x_{\beta }\right)\frac{\mathrm{d}{L}_{\beta i}\left(x_{\beta }\right)}{\mathrm{d}x}{\underset{¯}{u}}_{\beta i}^{\mathrm{ph}})}\\ =1·\sum _{i\in \ell \left(\beta \right)}\left(\frac{\mathrm{d}{L}_{\beta i}\left(x_{\beta }\right)}{\mathrm{d}x}{\underset{¯}{u}}_{\beta i}^0\right)\hfill \\ \multicolumn{1}{c}{=[\frac{{\underset{¯}{u}}_{\beta 2}^0}{h_{\beta }}+\frac{2{\underset{¯}{u}}_{\beta 3}^0}{h_{\beta }}·\left(\frac{x_{\beta }-x_{\beta }}{h_{\beta }}\right)+\frac{3{\underset{¯}{u}}_{\beta 4}^0}{h_{\beta }}·{\left(\frac{x_{\beta }-x_{\beta }}{h_{\beta }}\right)}^2]}\\ =\frac{{\underset{¯}{u}}_{\beta 2}^0}{h_{\beta }}\hfill \end{array}$$

(15)

Therefore, if $h_{\beta }$ is set to be 1, ${\underset{¯}{u}}_{\beta 2}^0$ has the physical meaning of the rotation at node β. Due to this feature, the proposed GFEM allows the straightforward imposition of displacement boundary conditions without using special techniques such as the penalty and Lagrange multiplier methods. ²¹

Numerical integration

As in other generalized/extended finite element applications with interfaces,^22,23 the numerical integration of elements with the plastic hinge enrichment defined by equation (8) requires the subdivision of integration elements. Figure 5(a) illustrates the numerical integration scheme for two-dimensional elements cut by a crack interface. ²³ In this case, the crack is described by Heaviside enrichment function, and the GFEM shape functions with this enrichment satisfy only C⁰ requirement at the location of crack interface. As a result, the standard integration rule such as Gauss quadrature cannot exactly integrate stiffness components computed with these shape functions if a single integration domain is used for the entire element. For the same reason, the fiber beam element with the plastic hinge enrichment needs to be subdivided at the plastic hinge location, at which only C¹ requirement is satisfied, and a numerical integration rule must be applied to each of the subdivided integration elements as illustrated in Figure 5(b). For the numerical integration of the GFEM discretization proposed in this study, the standard Gauss quadrature rule is adopted.

Figure 5.

Numerical integration scheme: (a) two-dimensional elements cut by a crack interface and (b) beam element with a plastic hinge enrichment.

Numerical experiments

In this section, the convergence analysis of the example shown in Figure 6 is performed to investigate the effectiveness of the proposed GFEM for the plastic hinge modeling based on the lumped plasticity. This model problem has a plastic hinge at the center of the span and is solved using both the GFEM and standard FEM for comparison purpose. The geometrical configuration of the model problem, material properties and magnitude of the linearly distributed loading are also indicated in the figure. The magnitude q₀ is 1.0.

Figure 6.

Model problem with an internal hinge at the center of the span subjected to linearly distributed loads.

The exact solution of this model problem (u_exact) is expressed by

$$u_{\mathrm{exact}}=\frac{q_0{x}^2}{120\mathrm{lEI}}(10l^3-10l^{2}x+5l{x}^2-x^3)(0\le x\le 2)$$

(16)

The exact strain energy of this problem can be calculated using the above exact solution, and is stated as

$$U_{\mathrm{exact}}=\frac{2{q_0}^2}{315l^2\mathrm{EI}}(35l^4-105l^3+126l^2-84l+24)$$

(17)

The relative error of convergence results $(e_r^E)$ is measured in energy norm and given in the form of

$$e_r^E=\sqrt{\frac{U_{\mathrm{exact}}-U_h}{U_{\mathrm{exact}}}}$$

(18)

where $U_h$ is the strain energy of the numerical solution obtained using the standard FEM or GFEM. This error estimate is used in the remaining part of this section.

Case with an even number of elements of equal size

In this section, a uniform h-extension is performed on the model problem using an even number of elements of equal size in the finite element mesh. As a result, the internal hinge is always located at the node of elements, as shown in Figure 7. In the standard FEM analysis, this can be modeled by changing the connectivity of the node represented by a circle in the figure.

Figure 7.

Uniform mesh refinement procedure for the standard FEM analysis with an even number of elements.

Table 1 presents the convergence analysis results obtained using the standard FEM, where the cubic Hermite shape functions (p = 3) are utilized. The number of elements, number of DOFs, relative error in energy norm, strain energy, and slope in log–log scale (ratio of the relative error in energy norm to the number of DOFs) are listed in the table. The results of the table are plotted in log–log scale in Figure 8 by having the number DOFs as x-axis and the relative error in energy norm as y-axis. As expected, the convergence rate of the standard FEM analysis approaches the theoretical value, ²¹ which is (p – 1) = 2, with increasing number of DOFs.

Table 1.

Convergence analysis results of h-extension obtained using the standard FEM with an even number of elements.

Number of elements	Numberof DOFs	Relative errorin energy norm	Strain energy	Slope in log–log scale
2	7	3.0000E–01	2.31111111111E–01	–
4	11	8.2916E–02	2.52222206208E–01	−2.845
8	19	2.1195E–02	2.53854164311E–01	−2.496
16	35	5.3274E–03	2.53961046047E-01	−2.260
32	67	1.3336E–03	2.53967802289E–01	−2.133
64	131	3.3353E–04	2.53968225716E–01	−2.067

DOF: degree of freedom.

Figure 8.

H-extension performed on the model problem using the standard FEM with an even number of elements.

Similarly, the results of uniform h-extension obtained using the GFEM with quadratic, cubic, and quartic polynomial enrichment functions are listed in Tables 2–4, respectively. The uniform mesh refinement procedure of this GFEM convergence analysis is illustrated in Figure 9, and the plastic hinge is modeled by the special enrichment functions given by equation (7), instead of the changing its nodal connectivity as in the standard FEM. The results of Tables 2–4 are plotted in log–log scale in Figure 10.

Table 2.

Convergence analysis results of h-extension obtained using the GFEM with an even number of elements and p = 2.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
2	12	0.0000E+00	2.53968253968E–01	–
4	18	2.0535E–01	2.43258762540E–01	N/A
8	30	2.4831E–01	2.38309115911E–01	0.372
16	54	1.4555E–01	2.48587986667E–01	−0.909
32	102	7.8247E–02	2.52413309712E–01	−0.976
64	198	4.0516E–02	2.53551353332E–01	−0.992

DOF: degree of freedom.

Table 3.

Convergence analysis results of h-extension obtained using the GFEM with an even number of elements and p = 3.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
2	16	0.0000E+00	2.53968253968E–01	–
4	24	0.0000E+00	2.53968253968E–01	N/A
8	40	2.0677E–02	2.53859672805E–01	N/A
16	72	6.0071E–03	2.53959089460E–01	−2.103
32	136	1.6034E–03	2.53967601043E–01	−2.077
64	264	4.1352E–04	2.53968210540E–01	−2.043

DOF: degree of freedom.

Table 4.

Convergence analysis results of h-extension obtained using the GFEM with an even number of elements and p = 4.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
2	20	0.0000E+00	2.53968253968E–01	–
4	30	0.0000E+00	2.53968253968E–01	N/A
8	50	7.7749E–04	2.53968100447E–01	N/A
16	90	1.1441E–04	2.53968250644E–01	−3.260
32	170	1.5308E–05	2.53968253909E–01	−3.163
64	330	2.1261E–06	2.53968253967E–01	−2.976

DOF: degree of freedom.

Figure 9.

Uniform mesh refinement procedure for the GFEM analysis with an even number of elements.

Figure 10.

H-extension performed on the model problem using the GFEM with an even number of elements.

Several interesting observations can be made from these results. First, linear convergence rates are achieved in all of the three cases. The theoretical convergence rate for uniform h-extension of smooth one-dimensional problems is given by (p – 1), where p is the polynomial order of shape functions. The convergence rates of the GFEM results with p = 2, 3, and 4 at the last convergence stage are 0.992, 2.043, and 2.976, respectively. This confirms that the convergence rates of the proposed GFEM analysis results coincide with the theoretical estimates very well, thus demonstrating its effectiveness and accuracy. Second, it can be seen from the results that the GFEM with quartic polynomial enrichments exhibits the lowest level of relative error for a given number of DOFs, thus confirming that it is the most effective one among the three cases discussed. Third, interestingly, the GFEM discretization can exactly approximate the solution of the problem if only two elements are used in the mesh. Similarly, the GFEM discretization with p = 3 or 4 can also produce the exact solution in case that four elements are used. The partition of unity functions in the proposed GFEM are cubic polynomials and thus the GFEM shape functions with high-order polynomial enrichments may be able to locally approximate the solution, of which polynomial order is even higher than that of the enrichment functions used, with high accuracy. Therefore, if only two or four elements are used in the domain, it may be possible to obtain a highly accurate solution using the GFEM shape functions. A similar observation was made by Son et al. ²⁴

Tables 5–7 summarize the results of p-extension performed using the GFEM shape functions with 4, 16, and 64 elements, respectively. The results of the tables are plotted in Figure 11 in log–log scale. These results indicate that exponential convergence rates can be achieved with increasing polynomial order of GFEM shape functions. This result coincides with theoretical estimates as discussed by Belytschko et al., ²² demonstrating the effectiveness and robustness of the proposed GFEM.

Table 5.

Convergence analysis results of p-extension obtained using the GFEM with four elements.

Polynomial order	Number of DOFs	Error in energy norm	Strain energy	Slope in log–log scale
Quadratic	18	2.0535E–01	2.43258762540E–01	–
Cubic	24	0.0000E+00	2.53968253968E–01	N/A
Quartic	30	0.0000E+00	2.53968253968E–01	N/A

DOF: degree of freedom.

Table 6.

Convergence analysis results of p-extension obtained using the GFEM with 16 elements.

Polynomial order	Number of DOFs	Error in energy norm	Strain energy	Slope in log–log scale
Quadratic	54	1.4555E–01	2.48587986667E–01	–
Cubic	72	6.0071E–03	2.53959089460E–01	−11.080
Quartic	90	1.1441E–04	2.53968250644E–01	−17.750

DOF: degree of freedom.

Table 7.

Convergence analysis results of p-extension obtained using the GFEM with 64 elements.

Polynomial order	Number of DOFs	Error in energy norm	Strain energy	Slope in log–log scale
quadratic	198	4.0516E–02	2.53551353332E–01	–
cubic	260	4.1352E–04	2.53968210540E–01	−16.830
quartic	325	2.1261E–06	2.53968253967E–01	−23.619

DOF: degree of freedom.

Figure 11.

P-extension performed on the model problem using the GFEM with an even number of elements.

Case with an odd number of elements of equal size

In this section, we perform convergence analyses for the model problem with an odd number of elements of equal size. As a result, the internal plastic hinge is always located inside an element, as shown in Figure 12. This case can be handled only using the proposed GFEM discretization, not by the standard FEM. In the GFEM discretization, the nodes of the element with the plastic hinge are enriched with the plastic hinge enrichment function described by equation (7).

Figure 12.

Uniform mesh refinement procedure for the GFEM analysis with an odd number of elements.

Tables 8–10 list the results of uniform h-extension obtained using the GFEM with quadratic, cubic, and quartic polynomial enrichment functions. The results of the tables are plotted in log–log scale in Figure 13. Similarly to the convergence analysis results of section “Computational cost analysis,” linear convergence rates are achieved in all of the three cases even though the internal plastic hinge is described by the plastic hinge enrichment function, not by an element node and its connectivity. The convergence rates of the GFEM results with p = 2, 3, and 4 at the last convergence stage are 1.013, 2.085, and 3.104, respectively. This again confirms that the convergence rates of the GFEM analysis results coincide with the theoretical estimates very well, thus demonstrating the effectiveness and accuracy of the proposed GFEM.

Table 8.

Convergence analysis results of h-extension obtained using the GFEM with an odd number of elements and p = 2.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
3	18	7.9584E–02	2.52359717319E–01	–
5	24	2.9694E–01	2.31575018768E–01	4.577
9	36	2.2904E–01	2.40645251657E–01	−0.640
17	60	1.3818E–01	2.49119057168E–01	−0.989
33	108	7.6038E–02	2.52499866046E–01	−1.016
65	204	3.9914E–02	2.53563650185E–01	−1.013

DOF: degree of freedom.

Table 9.

Convergence analysis results of h-extension obtained using the GFEM with an odd number of elements and p = 3.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
3	24	7.5664E–04	2.53968108570E–01	–
5	32	8.9971E–03	2.53947695795E–01	8.606
9	48	1.6632E–02	2.53898000400E–01	1.515
17	80	5.3531E–03	2.53960976335E–01	−2.219
33	144	1.5102E–03	2.53967674742E–01	−2.153
65	272	4.0094E–04	2.53968213142E–01	−2.085

DOF: degree of freedom.

Table 10.

Convergence analysis results of h-extension obtained using the GFEM with an odd number of elements and p = 4.

Number of elements	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
3	30	0.0000E + 00	2.53968253968E–01	–
5	40	3.6338E–05	2.53968253633E–01	N/A
9	60	5.1511E–04	2.53968186581E–01	6.539
17	100	9.2535E–05	2.53968251794E–01	−3.361
33	180	1.3602E–05	2.53968253921E–01	−3.262
65	340	1.8891E–06	2.53968253967E–01	−3.104

DOF: degree of freedom.

Figure 13.

H-extension performed on the model problem using the GFEM with an odd number of elements.

Tables 11–13 summarize the results of p-extension performed on the model problem using the GFEM shape functions with 5, 17, and 65 elements, respectively. The results of the tables are plotted in log–log scale in Figure 14. These results indicate that exponential convergence rates can be achieved with increasing polynomial order of GFEM shape functions. This result coincides with theoretical estimates as discussed by Szabo and Babuška, ²¹ demonstrating the effectiveness and robustness of the proposed GFEM even in the case where the internal plastic hinge is described by the special plastic hinge enrichment function.

Table 11.

Convergence analysis results of p-extension obtained using the GFEM with five elements.

P order	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
Quadratic	24	2.9694E–01	2.31575018768E–01	–
Cubic	32	8.9971E–03	2.49119057168E–01	−12.154
Quartic	40	3.6338E–05	2.53563650185E–01	−24.701

DOF: degree of freedom.

Table 12.

Convergence analysis results of p-extension obtained using the GFEM with 17 elements.

P order	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
Quadratic	60	1.3818E–01	2.53947695795E–01	–
Cubic	80	5.3531E–03	2.53960976335E–01	−11.300
Quartic	100	9.2535E–05	2.53968213142E–01	−18.185

DOF: degree of freedom.

Table 13.

Convergence analysis results of p-extension obtained using the GFEM with 65 elements.

P order	Number of DOFs	Relative error in energy norm	Strain energy	Slope in log–log scale
quadratic	204	3.9914E–02	2.53968253633E–01	–
cubic	272	4.0094E–04	2.53968251794E–01	−15.992
quartic	340	1.8891E–06	2.53968253967E–01	−24.010

DOF: degree of freedom.

Figure 14.

P-extension performed on the model problem using the GFEM with an odd number of elements.

Efficient first-order plastic hinge analysis based on the GFEM

In this section, the first-order plastic hinge analysis based on the proposed GFEM is discussed, and its efficient solution technique based on the static condensation scheme is proposed. The effectiveness of this approach is investigated by analyzing the results of a model problem solved using both of the standard FEM and proposed GFEM.

First-order plastic hinge analysis

The first-order plastic hinge analysis is a series of linear analyses, where a structural collapse is simulated by the successive formation of plastic hinges in frame members. ¹² This is one of the simplest types of analyses to simulate the collapse of a building structure, as it does not require the nonlinear iteration procedure, which is generally the most time-consuming part of nonlinear analyses.

In the first-order plastic hinge analysis, if internal bending moment reaches the limit capacity at a certain location, a plastic hinge is formed and its bending resistance is kept constant as the limit value while other plastic hinges develop, as shown in Figure 15. The plastic hinge modeling is based on the lumped plasticity and thus the plastic hinge can be formed by inserting an internal hinge (or pin) at its location. Then, this procedure is repeated until the entire structure becomes unstable. The typical procedure of the first-order plastic hinge analysis is summarized in Figure 16.

Figure 15.

Plastic hinge moment–rotation relation.

Figure 16.

Typical procedure of the first-order plastic hinge analysis.

Figure 17 shows an example of the first-order plastic hinge analysis performed on a single-span frame structure using the standard FEM. A linear elastic analysis is performed on this structure in order to find the location of the first plastic hinge, which is formed inside the horizontal frame member (member 3) in this example, as shown in Figure 17(c). In order to describe it, the member 3 needs to be replaced by two frame elements (members 5 and 6), and an internal hinge is inserted at the node between the two new elements (Figure 17(d)). Then, the second linear elastic analysis is performed on the modified structure by applying additional vertical and horizontal loads, and the second plastic hinge is formed at the right end of member 6 (Figure 17(e)). This plastic hinge is located at an elemental node, so it can be modeled by inserting an internal hinge here without adding more elements. With the addition of the second plastic hinge, the structure becomes unstable, and thus the progressive collapse analysis is over.

Figure 17.

First-order plastic hinge analysis using the standard FEM: (a) original structure, (b) deformation due to applied forces, (c) formation of the first plastic hinge, (d) insertion of a node at the location of the first plastic hinge, (e) formation of the second plastic hinge, and (f) connectivity change of the node at the location of the second plastic hinge.

From this description of the first-order plastic hinge analysis using the standard FEM, it is clear that the structure has to be modified after the formation of a plastic hinge and reanalyzed from scratch to find the next plastic hinge location, since it requires either the replacement/addition of elements or change in connectivity between elements. In general, a real-world example of frame structural collapses requires the formation of thousands of plastic hinges, and thus its computational cost is enormously high. Consequently, we may need a more efficient methodology that requires less computational cost.

Figure 18 illustrates the first-order plastic hinge analysis performed on the same example as in Figure 17 using the proposed GFEM. The main difference between the two approaches is the description of plastic hinges. In order to describe the first plastic hinge, the proposed GFEM enriches the two nodes of the horizontal element with the special plastic hinge enrichment discussed in section “First-order plastic hinge analysis,” instead of modifying the mesh of the frame structure as in the case of the standard FEM (Figure 18(d)). Similarly, the second plastic hinge can be inserted by enriching only the node at the second plastic hinge location with the plastic hinge enrichment (Figure 18(f)). Apparently, the plastic hinge enrichment of the proposed GFEM requires the addition of more DOFs to the existing GFEM discretization but does not require the modification of the existing FE mesh. This advantage can be greatly amplified with the aid of the solution technique based on the static condensation, which is introduced in detail in the next section.

Figure 18.

First-order plastic hinge analysis using the proposed GFEM: (a) original structure, (b) deformation due to applied forces, (c) formation of the first plastic hinge, (d) enrichment of the nodes of the element with the first plastic hinge, (e) formation of the first plastic hinge, and (f) enrichment of the node with the second plastic hinge.

Generalized finite element solution technique based on the static condensation

As already discussed in section “Generalized finite element discretization for plastic hinge modeling based on lumped plasticity,” any GFEM shape functions can be hierarchically added to an existing GFEM discretization. This guarantees that in the first-order plastic hinge analysis, the stiffness matrix of the problem at the previous step is always nested in that at the next step because the plastic hinges once added to the structure never disappear. This unique feature enables the application of the static condensation scheme proposed by Gupta et al. ¹³ in the first-order plastic hinge analysis performed using the GFEM discretization discussed in section “Plastic hinge modeling using the GFEM.”

If a plastic hinge is inserted in the frame structure using the proposed GFEM shape functions, the GFEM approximation can be expressed by equation (11). As already discussed in section “Generalized finite element discretization for plastic hinge modeling based on lumped plasticity,” this composition of the GFEM approximation can result in the stiffness matrix form given by

$$\underset{{\mathbf{K}}_G}{\underbrace{\left[\begin{array}{cc}{\mathbf{K}}^0\hfill & {\mathbf{K}}^{0,\mathrm{ph}}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}^{\mathrm{ph},0}}& {\mathbf{K}}^{\mathrm{ph}}\hfill \end{array}\right]}}\underset{{\underset{¯}{\mathbf{u}}}_G}{\underbrace{\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]}}=\underset{{\mathbf{F}}_G}{\underbrace{\left[\begin{array}{c}{\mathbf{F}}^0\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}^{\mathrm{ph}}}\end{array}\right]}}$$

(19)

In most of the building structures that are designed and constructed in industry, plastic hinges are generally formed relatively in a small number of frame elements at the final stage of structural failure in comparison with the total number of elements, as can be seen in the example of Figure 2. Consequently, in the first-order plastic hinge analysis based on the proposed GFEM, the number of DOFs associated with plastic hinge enrichments (p) is much smaller than that associated with polynomial enrichments (n). This condition can be stated as

$$\dim \left({\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}\right)<<\dim \left({\underset{¯}{\tilde{\mathbf{u}}}}^0\right)$$

(20)

Therefore, if the factorized form of ${\mathbf{K}}^0$ is preserved at the initial step of the first-order plastic hinge analysis, the solutions at subsequent steps can be obtained only with a small additional computational cost by utilizing it. For this, we apply the static condensation scheme introduced by Jukić et al. ⁸ The term ${\underset{¯}{\tilde{\mathbf{u}}}}^0$ of equation (19) can be written as

$${\underset{¯}{\tilde{\mathbf{u}}}}^0={\left({\mathbf{K}}^0\right)}^{-1}{\mathbf{F}}^0-{\left({\mathbf{K}}^0\right)}^{-1}{\mathbf{K}}^{0,\mathrm{ph}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}={\underset{¯}{\mathbf{u}}}^0-{\mathbf{S}}^{0,\mathrm{ph}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$$

(21)

where ${\mathbf{S}}^{0,\mathrm{ph}}$ is the pseudo solutions, which is defined by

$${\mathbf{S}}^{0,\mathrm{ph}}:={\left({\mathbf{K}}^0\right)}^{-1}{\mathbf{K}}^{0,\mathrm{ph}}$$

(22)

If the factorized form of ${\mathbf{K}}^0$ is preserved at the initial step of the analysis, the pseudo solutions can be computed through forward- and back-substitutions on ${\mathbf{K}}^0$ as shown below

$$\underset{\mathrm{pseudo}\mathrm{solutions}}{\underbrace{{{\mathbf{K}}^0\mathbf{S}}^{0,\mathrm{ph}}}}=\underset{\mathrm{pseudo}\mathrm{loads}}{\underbrace{{\mathbf{K}}^{0,\mathrm{ph}}}}$$

(23)

Then, the terms ${\mathbf{K}}^{\mathrm{ph}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ of equation (19) can be stated as

$${\mathbf{K}}^{\mathrm{ph}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}={\mathbf{F}}^{\mathrm{ph}}-{\mathbf{K}}^{\mathrm{ph},0}[{\underset{¯}{\mathbf{u}}}^0-{\mathbf{S}}^{0,\mathrm{ph}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}]$$

(24)

This can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{\underset{{\hat{\mathbf{K}}}^{\mathrm{ph}}}{\underbrace{[{\mathbf{K}}^{\mathrm{ph}}-{\mathbf{K}}^{\mathrm{ph},0}{\mathbf{S}}^{0,\mathrm{ph}}]}}{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}=\underset{{\hat{\mathbf{F}}}^{\mathrm{ph}}}{\underbrace{{\mathbf{F}}^{\mathrm{ph}}-{\mathbf{K}}^{\mathrm{ph},0}{\underset{¯}{\mathbf{u}}}^0}}}\end{array}$$

(25)

The dimensions of matrices ${\hat{\mathbf{K}}}^{\mathrm{ph}}$ , ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ , and ${\hat{\mathbf{F}}}^{\mathrm{ph}}$ are(p × p), (p × 1), and (p × 1), respectively. The size of the stiffness matrix system given by equation (25) is smaller than that of the original system given by equation (19), its solution ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ can be directly obtained only with a very small computational cost. By plugging the obtained ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ into equation (21), we can easily compute ${\underset{¯}{\tilde{\mathbf{u}}}}^0$ . Then, the solution at any subsequent step with plastic hinges formed in the structure can be constructed from

$$u_h={\tilde{u}}^0+u^{\mathrm{ph}}=\left[\begin{array}{cc}{\mathbf{N}}^0\hfill & {\mathbf{N}}^{\mathrm{ph}}\hfill \end{array}\right]\left[\begin{array}{c}{\underset{¯}{\tilde{\mathbf{u}}}}^0\hfill \\ \multicolumn{1}{c}{{\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}}\end{array}\right]$$

(26)

In the procedure described up to this point, most of the computational cost is spent in the computations stated by equations (22) and (25). This is further discussed in the next section.

Numerical experiments

This section investigates the effectiveness of the GFEM solution technique discussed in section “Generalized finite element solution technique based on the static condensation” by solving a model problem illustrated in Figure 19 using both the standard FEM and proposed GFEM and by comparing their computational costs. The model problem shown in the figure is basically a multi-span continuous beam, of which spans have a constant length (L) and bending stiffness (EI). Its spans can be categorized into two types such as those with and without concentrated loads. The spans with a single concentrated load (P) represent the frame members with plastic hinges, while the free spans represent those without plastic hinges. The number of the concentrated loads in the continuous beam is denoted by n_A, and the number of the free spans is denoted by n_B. The limit bending capacity for plastic hinge formation is set to be the plastic bending moment. All of the blue and orange elements in Figure 19 have the same plastic bending moments, respectively, and the plastic bending moment of the former $(M_{p,\mathrm{blue}})$ is 10 times higher than that of the latter $(M_{p,\mathrm{orange}})$ . The following parameters are adopted for the analysis: P = –1.0, EI = 100.0, L = 4.0, $M_{p,\mathrm{orange}}=10.0$ .

Figure 19.

Model problem for the cost analysis of the proposed GFEM.

In order to avoid the simultaneous formation of many plastic hinges, the frame members with orange color have ten times greater plastic moment than those with blue color. Under this condition, (2n_A – 1) plastic hinges are formed in the structure before it becomes completely unstable. The progressive collapse procedure of this structure with n_A = 5 is illustrated in Figure 20, and the plastic hinges are represented by red solid circles. In this case, nine plastic hinges are formed in the structure before its collapse, as indicated in the figure.

Figure 20.

Progressive collapse mechanisms of the model problem.

Computational cost analysis

In this section, the computational costs required to perform the first-order plastic hinge analysis of the model problem using the standard FEM and proposed GFEM are defined in terms of an accumulated CPU time, respectively, and their computational costs are compared. If its first-order plastic hinge analysis is performed using the standard FEM, a direct solver such as Cholesky factorization needs to be used at all steps. As a result, its accumulated CPU time required to perform the analysis up to step (2n_A – 1) (T_FEM) can be estimated by the following equation

$$T_{\mathrm{FEM}}=\sum _{i=1}^{(2N_A-1)}{t}_{\mathrm{FEM}}^i$$

(27)

where $t_{\mathrm{FEM}}^i$ is the CPU time required to compute the solution of the linear system of equations constructed with the standard FEM approximation at the first step i using the direct solver.

In the proposed GFEM plastic hinge analysis, the direct solver is used only at the first step, and the solutions at subsequent steps can be obtained using the static condensation scheme discussed in section “Generalized finite element solution technique based on the static condensation.” In the procedure of the static condensation, most of the computational cost is spent during the computations expressed by equations (22) and (25), as already mentioned in section “Generalized finite element solution technique based on the static condensation.” Therefore, the accumulated CPU time required by the GFEM analysis up to step (2n_A – 1) (T_GFEM) can be estimated as

$$T_{\mathrm{GFEM}}=t_{\mathrm{GFEM}}^1+\sum _{i=2}^{(2N_A-1)}(t_S^i+t_{\mathrm{ph}}^i)$$

(28)

where $t_{\mathrm{GFEM}}^1$ is the CPU time required to compute the solution of the linear system of equations built with the proposed GFEM approximation at the first step using the direct solver. In the above equation, $t_S^i$ and $t_{\mathrm{ph}}^i$ are the CPU times required to construct ${\mathbf{S}}^{0,\mathrm{ph}}$ in equation (22) and to compute ${\underset{¯}{\mathbf{u}}}^{\mathrm{ph}}$ of equation (25) at step i, respectively.

The first-order plastic hinge analyses of the model problem are performed under three different conditions such as n_B = 800, n_B = 1600 and n_B = 3200 using both the standard FEM and proposed GFEM with cubic polynomial enrichments. Here, in all of the three cases, n_A is constantly equal to 5, as in the case shown in Figure 19. The number of DOFs of the standard FEM and proposed GFEM approximations at the first step are summarized in Table 14. This table indicates that the GFEM approximation requires twice more the DOFs than the standard FEM approximation. As can be seen from equation (4), the cubic GFEM shape functions require four DOFs per node, while the cubic Hermite shape functions of the standard FEM requires only two.

Table 14.

Number of DOFs required by the standard FEM and GFEM approximations at the first step.

nB	Standard FEM	GFEM
800	1630	3260
1600	3230	6460
3200	6430	12,860

FEM: finite element method; GFEM: generalized finite element method.

The computational costs of the standard FEM and proposed GFEM analyses for the three cases are plotted in Figures 21–23, respectively. Several interesting observations can be made from these results. First, in the standard FEM analysis, the accumulated CPU time increases in proportion with the number of plastic hinges formed in the structure. In contrast, the GFEM analysis requires a relatively high amount of CPU time to obtain the first step solution, but from the second step, the required CPU time becomes much less and gradually increases as the number of plastic hinges increases. This confirms that the static condensation scheme proposed in section “Generalized finite element solution technique based on the static condensation” can greatly reduce the computational cost to obtain the solutions after the first step. As a consequence, the accumulated CPU time required for the GFEM analysis becomes smaller than that of the standard FEM analysis if a sufficient number of plastic hinges are formed in the structure before its collapse.

Figure 21.

Computational costs of the standard FEM and GFEM analyses for n_B = 800.

Figure 22.

Computational costs of the standard FEM and GFEM analyses for n_B = 1600.

Figure 23.

Computational costs of the standard FEM and GFEM analyses for n_B = 3200.

Second, the cost effectiveness of the proposed GFEM becomes better if the size of the problem increases. For example, it can be seen from Figure 21 that the accumulated CPU time required by the GFEM analysis is larger than that by the standard FEM analysis, but the computational cost of the GFEM analysis becomes smaller than that of the standard FEM when the number of plastic hinges is more than eight in Figures 22 and 23.

Third, in the number of plastic hinges versus the accumulated CPU time graph shown in Figures 21–23, the slope of the GFEM analysis gradually increases as the number of steps increases. This happens mainly because the static condensation scheme proposed in section “Generalized finite element solution technique based on the static condensation” may not be computationally effective if the condition stated by equation (20) is not fully satisfied. In other words, it becomes less effective with increasing number of plastic hinges. In order to address this issue, the alternating scheme of direct solution and static condensation is proposed in the next section.

Alternating scheme of direct solution and static condensation

Figure 24 illustrates the alternating scheme of direct solution and static condensation in order to alleviate the effectiveness degradation of the static condensation scheme with increasing number of plastic hinges formed in the structure. It can be seen from the comparison of equations (27) and (28) that the static condensation scheme may require more CPU time than the direct solver if $(t_S^i+t_{\mathrm{ph}}^i)$ becomes greater than $t_{\mathrm{FEM}}^i$ . Therefore, we can solve the problem using the direct solver from scratch if this condition is met, and start the static condensation scheme again by utilizing the factorized form of the stiffness matrix, which is kept in the direction solution procedure at the previous step. In fact, it is difficult to obtain $t_{\mathrm{FEM}}^i$ unless the problem is solved using the standard FEM. Alternatively, $t_{\mathrm{GFEM}}^i/2^3$ may be used, since the number of DOFs required by the GFEM approximation is twice larger than that of the standard FEM, and the computational cost of a direct solver is generally proportional to the cube of the number of DOFs.

Figure 24.

Alternating scheme of direct solution and static condensation.

As an example to demonstrate the effectiveness of the proposed alternating scheme, the model problem of a larger size is solved for n_A = 50 and n_B = 3200 using both the standard FEM and GFEM with the alternating scheme. In this case, a total of 99 plastic hinges are formed before the collapse of the structure. The number of DOFs at the first step required by the standard FEM and GFEM are 6550 and 13,100, respectively. The computational costs of the two methods are plotted in Figure 25. It can be noted from the figure that the accumulated CPU time of the GFEM approach is greatly reduced with the aid of the alternating scheme, and the difference between the costs of two methods becomes more significant as the number of plastic hinges formed in the structure increases. This confirms that the proposed alternating scheme is highly effective to enhance the solution technique based on static condensation.

Figure 25.

Cost reduction of the GFEM analyses performed with the proposed alternating scheme for n_A = 50 and n_B = 3200.

Conclusion

In this study, we proposed a generalized finite element formulation for plastic hinge modeling based on lumped plasticity in the classical Euler–Bernoulli beam. In this approach, the plastic hinges are modeled using a special enrichment function, which can describe the weak discontinuity of the solution at the location of the plastic hinge. Furthermore, it is also possible to insert a plastic hinge at an arbitrary location of the element without modifying its connectivity or adding more elements. Instead, the formations of the plastic hinges are achieved by hierarchically adding more DOFs to existing elements. Due to these features, the proposed methodology can efficiently perform the first-order plastic hinge analysis of large-frame structures. A generalized finite element solution technique based on the static condensation scheme was also proposed in order to reduce the computational cost of a series of linear elastic problems, which is in general the most time-consuming portion of the first-order plastic hinge analysis. The main conclusions of this study can be summarized as follows:

1. The proposed GFEM is able to achieve the theoretical convergence rates of the static analysis, which are in principle identical to those of the standard FEM, for various polynomial orders of its shape functions such as quadratic, cubic, and quartic orders in both cases, where plastic hinges are located at nodes and within an element.

2. The GFEM approximations with high-order polynomial enrichments such as p = 4 are effective in achieving the low level of relative error in energy norm for a given number of DOFs.

3. The GFEM shape functions with high-order polynomial enrichments are able to obtain an almost exact solution if only one or two elements are used in the domain. This seems possible because the partition of unity functions used in this approach are the cubic Hermite shape functions associated with displacement, and the GFEM shape functions with high-order polynomial enrichments may be able to locally approximate the solution, of which polynomial order is even higher than that of the enrichment functions used.

4. The solution technique based on the static condensation is able to significantly reduce the computational cost of the first-order plastic hinge analysis using the GFEM approximation proposed in this study. Its effectiveness is especially high when the size of the problem is very large than the total number of plastic hinges formed in the structure before its collapse.

5. The alternating scheme of direct solution and static condensation is able to prevent the effectiveness degradation of the static condensation that may occur if the number of plastic hinges formed in the structure increases.