HDQ-based solution strategy for in-plane nonlinear buckling of shallow steel arches with linear elastic supports

This paper investigates the in-plane nonlinear buckling of shallow steel arches with linear elastic supports. Differential equations of equilibrium are built considering a virtual work principle. To balance computational accuracy, efficiency, and cost, a harmonic differential quadrature method (HDQ)-based solution strategy combined with an iterative process is utilized to follow the complete buckling path. The corresponding critical buckling load is obtained. The convergence of HDQ algorithm is then examined, and the accuracy of presented solutions is validated by comparing the FE predictions. It was found that the nonlinear behavior of circular arches with linear elastic supports can be precisely predicted by the presented HDQ-based solution strategy. The arch could buckle either in a limit point mode or a bifurcation mode, depending on the modified slenderness ratio.

1 Introduction

Arches are commonly used in civil, mechanical, and aerospace engineering [1]. It can convert external force into inner axial force and bending moment. Due to its unique bearing characteristics, arches may collapse abruptly within elastic range, leading to structural failure [2]. This complex mechanical behavior has attracted widespread attention from scholars both domestically and internationally.

Considering that ideal fixed and hinged boundaries do not exist in practical engineering, the actual boundary is usually approximated by linear springs, torsional springs, elastic foundations, etc. In accordance with the principle of stiffness equivalence. Pi and Bradford [3, 4] investigated the nonlinear in-plane buckling of circular arches with linear elastic supports and rotational end restraints. Closed form solutions for the critical buckling load were derived by using an energy principle. Internal connections between restraint stiffness and buckling shape were revealed. Han et al. [5] proposed a 9-node assumed natural strain shell element dealing with the geometrical nonlinear buckling of deep circular arches with arbitrary elastic edge supports. The primary path at any point on the axis during buckling were recorded. The in-plane snap-through buckling of sinusoidal arches with rotational supports that stiffen under compressed was analyzed by plaut [6]. He claimed that the critical buckling load of arch may be increased significantly due to support stiffening, especially if the initial flexural resistance is small. Cai and his-corkers [7, 8] numerically explored the snap-through behavior of elastic shallow inextensible circular arches with variable elastic horizontal supports under unilateral displacement control. They found the critical stiffness increases with the decrease of arch length, and the decrease of horizontal stiffness can expand the snap region. Zhou et al. [9] confirmed that shallow arches with symmetric elastic supports can bifurcate into secondary paths with high-order symmetric modes under the action of uniformly distributed vertical load, from a theoretical perspective. Han et al. [10, 11] carried out elastic-plastic in-plane buckling experiments of steel circular arch with horizontal elastic supports and torsional constraints. The complete equilibrium path, as well as the corresponding failure mode were captured. Design method for the ultimate bearing capacity of arches has been proposed for engineering applications. Hu et al. [12, 13] explored the nonlinear elastic stability of pin-ended parabolic multi-span continuous arches. In their analysis, the horizontal connections between loaded and unloaded elements were simulated by linear springs. Despite extensive research have been conducted and innovative results have been found, there are still many inconveniences in determining the critical buckling load of arches under complex boundary conditions in open literature, such as complex formulas and high computational costs.

To fulfill the research gap, this paper proposes a HDQ-based solution strategy dealing with the in-plane nonlinear buckling of shallow steel arches with linear elastic support. The differential equations of equilibrium were built initially utilizing a virtual work principle. By introducing the differential quadrature element (HDQ) algorithm, the governing equations were discretized. An iterative process was then performed to follow the limit point buckling path. After that, the displacement perturbation was conducted to establish the governing equations for bifurcation, the corresponding path was traced. Finally, the convergence of HDQ algorithm and the accuracy of the presented solutions were examined.

2 Limit point buckling analysis

Consider a steel arch having a bi-symmetric I-shape cross-section subjected to uniformly distributed radial loads, as shown in Figure 1. The Young’s modulus of the steel $E=2\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$, Poisson’s ratio $\mu =0.3$. The radius and included angle of the arch are denoted as $R$ and $2\mathrm{\Theta }$. The overall height (width) of the cross-section $h=400\text{\hspace{0.17em}mm}$ ($b=180\text{\hspace{0.17em}mm}$), and the flange (web) thickness $t_f=13.5\text{\hspace{0.17em}mm}$ ($t_w=8.6\text{\hspace{0.17em}mm}$). The $x$ and $y$ axes are arranged at the principal centroid axes of the cross-section.

Figure 1

Figure 1. Geometry model of shallow steel arches under uniformly distributed radial load.

It is assumed the circular arch satisfies the classic Euler Bernoulli deformation condition. In the framework of virtual work principle, the variation of the total potential energy can be expressed as [14].

$\begin{array}{ll}\hfill & \hfill \delta \mathrm{\Pi }={\int }_0^{2\mathrm{\Theta }}\left\{NR\delta \left[{\stackrel{\sim }{w}}^{\prime }-\stackrel{\sim }{v}+\frac{1}{2}{\left({\stackrel{\sim }{v}}^{\prime }\right)}^2\right]+M\delta {\stackrel{\sim }{v}}^{″}+q{R}^2\delta \stackrel{\sim }{v}\right\}\mathrm{d}\theta \\ \hfill & +{\displaystyle \sum _{i=1}^2}\left(k_{v_i}\stackrel{\sim }{v}\delta {\stackrel{\sim }{v}}_i+k_{w_i}\stackrel{\sim }{w}\delta {\stackrel{\sim }{w}}_i\right)R^2=0\hfill \end{array}(1)$

where $\stackrel{\sim }{v}=v/R$ and $\stackrel{\sim }{w}=w/R$ are the dimensionless radial and axial displacements, $R$ is the radius of the arch, $q$ is the load density, ${\left(\right)}^{\prime }=\mathrm{d}\left(\right)/\mathrm{d}\theta$, ${\left(\right)}^{″}={\mathrm{d}}^2\left(\right)/\mathrm{d}{\theta }^2$, and $\theta$ is the angular coordinate, $k_v$ and $k_w$ are the axial and radial stiffness coefficients of the elastic supports, and number 1 (2) at the upper (lower) limit of cumulative symbol represents the left (right) arch foot, $N$ and $M$ are the axial compressive forces and bending moment, which are defined by

$N=-{\int }_A\sigma \mathrm{d}A=-EA\left[{\stackrel{\sim }{w}}^{\prime }-\stackrel{\sim }{v}+\frac{1}{2}{\left({\stackrel{\sim }{v}}^{\prime }\right)}^2\right](2)$

$M={\int }_A\sigma y\mathrm{d}A=-\frac{EI}{R}{\stackrel{\sim }{v}}^{″}(3)$

where $A$ and $I$ are the cross-sectional area and second moment of area of cross-section about the major principal axis, $\sigma$ denotes the cross-section normal stress and $E$ is the Young’s modulus.

Integrating Equation 1 by parts leads to [15].

$\delta \stackrel{\sim }{w}: N^{\prime }=0\Rightarrow N=\text{const}(4)$

$\delta \stackrel{\sim }{v}: M^{″}+NR{\stackrel{\sim }{v}}^{″}+NR=q{R}^2(5)$

Substituting Equations 2, 3 into Equations 4, 5 has

${\stackrel{\sim }{w}}^{″}-{\stackrel{\sim }{v}}^{\prime }+{\stackrel{\sim }{v}}^{\prime }{\stackrel{\sim }{v}}^{″}=0(6)$

$\frac{EI}{R}{\stackrel{\sim }{v}}^{\left(iv\right)}+EAR\left({\stackrel{\sim }{w}}^{″}-{\stackrel{\sim }{v}}^{\prime }+{\stackrel{\sim }{v}}^{\prime }{\stackrel{\sim }{v}}^{″}\right)+EAR\left[{\stackrel{\sim }{w}}^{\prime }-\stackrel{\sim }{v}+\frac{1}{2}{\left({\stackrel{\sim }{v}}^{\prime }\right)}^2\right]{\stackrel{\sim }{v}}^{″}=-q{R}^2(7)$

To obtain the real root of Equations 6, 7 a harmonic differential quadrature algorithm is adopted to solve the governing equation in spatial domain. The $\xi \text{th}$-order derivation of the function $\stackrel{\sim }{w}\left(\stackrel{\sim }{v}\right)$ at any discrete point can be expressed as [16]

$\left\{\stackrel{\sim }{v},\stackrel{\sim }{w}\right\}={\displaystyle \sum _{j=1}^n}{h}_j\left(\xi \right)\left\{{\stackrel{\sim }{v}}_m,{\stackrel{\sim }{w}}_m\right\}\text{\hspace{0.5em}and\hspace{0.17em}}\frac{{\partial }^m}{\partial {\xi }^m}{\left\{\stackrel{\sim }{v},\stackrel{\sim }{w}\right\}}_{\xi ={\xi }_j}={\displaystyle \sum _{j=1}^n}{C}_{ij}^{\left(m\right)}\left\{{\stackrel{\sim }{v}}_m,{\stackrel{\sim }{w}}_m\right\}(8)$

where $m$ is the number of sampling points, $h_i\left(\xi \right)$ is the harmonic test function, and is defined by [17]

$h_i\left(\xi \right)={\displaystyle \prod _{k=0,k\ne i}^n}\sin \left[\pi \left(\xi -{\xi }_k\right)/2\right]/{\displaystyle \prod _{k=0,k\ne i}^n}\sin \left[\pi \left({\xi }_i-{\xi }_k\right)/2\right](9)$

The first and second derivatives of weighting coefficients for $i\ne j$ can be written using the following formula [18]:

$C_{ij}^{\left(1\right)}=\frac{\pi P\left({\xi }_i\right)/2}{P\left({\xi }_j\right)\sin \left[\pi \left({\xi }_i-{\xi }_j\right)/2\right]}i,j=1,2\dots n(10)$

$C_{ij}^{\left(2\right)}=C_{ij}^{\left(1\right)}\left\{2C_{ij}^{\left(1\right)}-\pi \cot \left[\frac{\pi \left({\xi }_i-{\xi }_k\right)}{2}\right]\right\}i,j=1,2\dots n(11)$

where

$P\left({\xi }_i\right)={\displaystyle \prod _{j=1,j\ne i}^n}\sin \left[\pi \left({\xi }_i-{\xi }_k\right)/2\right]i,j=1,2\dots n(12)$

The weighting coefficients of the first-order and second-order derivatives for $i=j$ can be written as

$C_{ii}^{\left(r\right)}=-{\displaystyle \sum _{j=1,j\ne i}^n}{C}_{ij}^{\left(r\right)}r=1\text{\hspace{0.17em}or\hspace{0.17em}}2\text{\hspace{0.17em}with\hspace{0.5em}}i=1,2\dots n(13)$

The weighting coefficient of the third- and fourth-order derivatives can be computed easily from $C_{ij}^{\left(1\right)}$ and $C_{ij}^{\left(2\right)}$ as

$C_{ij}^{\left(3\right)}={\displaystyle \sum _{k=1}^n}{C}_{ij}^{\left(1\right)}{C}_{ij}^{\left(2\right)}\text{\hspace{0.5em}and\hspace{1em}}{C}_{ij}^{\left(4\right)}={\displaystyle \sum _{k=1}^n}{C}_{ij}^{\left(2\right)}{C}_{ij}^{\left(2\right)}(14)$

The Chebyshev–Gauss–Lobatto polynomial with adjacent-δ points is used to generate the nodal grid point [19]:

$\begin{array}{ll}\hfill & \hfill {\xi }_1=0,{\xi }_2={10}^{-4},{\xi }_i\\ \hfill & =\frac{1}{2}{\left\{1-\cos \left[\frac{\left(i-1\right)\pi }{N-1}\right]\right\}}_{i=3,4,\dots k-2},{\xi }_{k-1}=1-{10}^{-4},{\xi }_k=1\hfill \end{array}(15)$

Considering Equations 8–15 and rewriting Equations 6, 7 in a differential quadrature form has

$\left\{{\mathit{K}}_L+{\mathit{K}}_{NL}\right\}\left\{\mathit{d}\right\}=-\left\{\mathit{Q}\right\}(16)$

with

$\left\{\mathit{d}\right\}={\left\{{\left\{{\stackrel{\sim }{w}}_i\right\}}^T,{\left\{{\stackrel{\sim }{v}}_i\right\}}^T\right\}}^T\text{\hspace{0.17em}and\hspace{0.17em}}\left\{\mathit{Q}\right\}={\left\{{\left\{\mathbf{0}\right\}}^T,{\left\{{qR}^2\right\}}^T\right\}}^T(17)$

${\mathit{K}}_L=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\Theta }}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}& -{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}\\ \frac{EAR}{4{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}& \frac{EI}{16R{\mathrm{\Theta }}^4}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(4\right)}-\frac{EAR}{2\mathrm{\Theta }}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}\end{array}\right\}(18)$

${\mathit{K}}_{NL}=\left\{\begin{array}{ll}\mathbf{0}& \frac{1}{4{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}\\ \frac{EAR}{32{\mathrm{\Theta }}^5}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}{\stackrel{\sim }{v}}_m& \begin{array}{c}\frac{EAR}{8{\mathrm{\Theta }}^3}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}-\frac{EAR}{16{\mathrm{\Theta }}^4}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}{\stackrel{\sim }{v}}_m\\ +\frac{EAR}{128{\mathrm{\Theta }}^6}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}\end{array}\end{array}\right\}(19)$

where ${\mathit{K}}_L$ is linear stiffness matrix, ${\mathit{K}}_{NL}$ is nonlinear stiffness matrix, $\left\{\mathit{Q}\right\}$ and $\left\{\mathit{d}\right\}$ are the displacement and load vectors.

The radial and axial elastic supports for circular arches can be written as

$EA\left[{\stackrel{\sim }{w}}_1-{\stackrel{\sim }{v}}_1+\frac{1}{8{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m\right]-k_{w_1}{\stackrel{\sim }{w}}_{1}R=0(20)$

$EA\left[{\stackrel{\sim }{w}}_n-{\stackrel{\sim }{v}}_n+\frac{1}{8{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(1\right)}{\stackrel{\sim }{v}}_m\right]+k_{w_n}{\stackrel{\sim }{w}}_{n}R=0(21)$

$\begin{array}{ll}\hfill & \hfill {\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(2\right)}{\stackrel{\sim }{v}}_m=\frac{1}{2\mathrm{\Theta }}EA\left[{\stackrel{\sim }{w}}_1-{\stackrel{\sim }{v}}_1+\frac{1}{8{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m\right]\\ \hfill & \times {\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m-\frac{EI}{{8R}^2{\mathrm{\Theta }}^3}{\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(3\right)}{\stackrel{\sim }{v}}_m+k_{v_1}{\stackrel{\sim }{v}}_{1}R=0\hfill \end{array}(22)$

$\begin{array}{ll}\hfill & \hfill {\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(2\right)}{\stackrel{\sim }{v}}_m=\frac{1}{2\mathrm{\Theta }}EA\left[{\stackrel{\sim }{w}}_n-{\stackrel{\sim }{v}}_n+\frac{1}{8{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(1\right)}{\stackrel{\sim }{v}}_m{\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(1\right)}{\stackrel{\sim }{v}}_m\right]\\ \hfill & \times {\displaystyle \sum _{m=1}^n}{C}_{1m}^{\left(1\right)}{\stackrel{\sim }{v}}_m-\frac{EI}{{8R}^2{\mathrm{\Theta }}^3}{\displaystyle \sum _{m=1}^n}{C}_{nm}^{\left(3\right)}{\stackrel{\sim }{v}}_m-k_{v_n}{\stackrel{\sim }{v}}_{n}R=0\hfill \end{array}(23)$

Since the nonlinear equations are hard to solve, an iterative method is employed to gradually approach the true solution. Ignoring the nonlinear stiffness matrix, Equation 16 can be further simplified as

$\left\{{\mathit{K}}_L\right\}\left\{\mathit{d}\right\}=-\left\{\mathit{Q}\right\}(24)$

Keeping Equations 17, 18, 20–24 in mind and performing

${\left\{\mathit{d}\right\}}_l\Rightarrow {\left\{{\mathit{K}}_L\right\}}^{-1}\left\{\mathit{Q}\right\}={\left\{{\left\{{\stackrel{\sim }{w}}_{il}\right\}}^T,{\left\{{\stackrel{\sim }{v}}_{il}\right\}}^T\right\}}^{T}i=1,2,\dots n(25)$

linear solution for the nodal deformation can be determined.

Substituting ${\left\{\mathit{d}\right\}}_l$ as an initial value into Equation 25 and solving Equations 16–19, the first-teration solution for the nonlinear displacement is obtained. Repeating iteration process until $\max \left\{\left.{\left\{\mathit{d}\right\}}_{k+1}-{\left\{\mathit{d}\right\}}_k\right\}\right\}<{10}^{-4}$, the nonlinear nodal displacement vector $\left\{{\mathit{d}}_k\right\}$ under $\left\{{\mathit{Q}}_{\mathit{k}}\right\}$ can be derived. Traversing the load density $q\in \left[0,1.5q_{cr}^l\right]$ subsequently, the complete limit point buckling path can be traced.

3 Bifurcation buckling analysis

Since the bifurcation point is in a primary equilibrium path $\left\{\stackrel{\sim }{w},\stackrel{\sim }{v}\right\}$ as well as in a bifurcation path $\left\{\stackrel{\sim }{w}+{\stackrel{\sim }{w}}_b,\stackrel{\sim }{v}+{\stackrel{\sim }{v}}_b\right\}$ with constant cross-sectional inner forces $\left\{N,M\right\}$ and external load $\left\{q\right\}$, the axial and radial governing equations for bifurcation can be stated as [20]

$N_b^{\prime }=0\Rightarrow {\stackrel{\sim }{w}}_b^{″}-{\stackrel{\sim }{v}}_b^{\prime }+{\stackrel{\sim }{v}}_b^{\prime }{\stackrel{\sim }{v}}_b^{″}=0(26)$

$EI{\stackrel{\sim }{v}}_b^{\left(iv\right)}+N_b{R}^2{\stackrel{\sim }{v}}_b^{″}=0(27)$

where subscript ${\left(\right)}_b=\delta \left(\right)$.

Writing Equations 26, 27 into a differential quadrature form has

$\left\{{\mathit{K}}_b\right\}\left\{{\mathit{d}}_{\mathit{b}}\right\}=\left\{\mathbf{0}\right\}(28)$

in which

$\left\{{\mathit{d}}_{\mathit{b}}\right\}={\left\{{\left\{{\stackrel{\sim }{w}}_{b_i}\right\}}^T,{\left\{{\stackrel{\sim }{v}}_{b_i}\right\}}^T\right\}}^T(29)$

${\mathit{K}}_b=\left\{\begin{array}{ll}\frac{1}{4{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}& \frac{1}{4{\mathrm{\Theta }}^3}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}{\stackrel{\sim }{v}}_{b_m}-\frac{1}{2\mathrm{\Theta }}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(1\right)}\\ \mathbf{0}& \frac{EI}{16{\mathrm{\Theta }}^4}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(4\right)}+\frac{N_b{R}^2}{4{\mathrm{\Theta }}^2}{\displaystyle \sum _{m=1}^n}{C}_{im}^{\left(2\right)}\end{array}\right\}(30)$

Substituting nonlinear nodal displacement vector $\left\{{\mathit{d}}_k\right\}$ and axial force $N_k$ into Equations 28–30 and checking whether $R\left({\mathit{K}}_b\right)<2n$ is satisfied, the corresponding bifurcation path can be followed.

Figure 2 illustrates the typical limit point path and bifurcation path of circular arches, according to Equations 16, 28. The horizontal axis takes the dimensionless radial displacement at crown, and vertical axis takes the dimensionless critical buckling load. Where $N_{TB}={\pi }^{2}EI/{\left(S/2\right)}^2$ represents the lowest load for limit point buckling, proposed by Pi and Bradford [1, 2]. $\lambda =\mathrm{\Theta }S/2r_x$ denotes the modified slenderness ratio, the radius of gyration $r_x=\sqrt{EI/EA}$, and ${\alpha }_v=4EI/k{S}^3$ (${\alpha }_w=EA/k{S}^3$) is the dimensionless flexible coefficient of the radial (axial) elastic supports. As observed, the buckling path exhibits highly nonlinear characteristics during the whole deformation processes. The bifurcation point may occur either before or after the upper limit point on the primary equilibrium path, depending on the value of the modified slenderness ratio.

Figure 2

Figure 2. In-plane buckling behaviors of circular arches: (a) limit point buckling path, (b) bifurcation buckling path.

4 Convergence and validation

To examine the convergence of HDQ algorithm in handling nonlinear problems, Table 1 lists the variations of critical buckling load of circular arches including limit point buckling as well as the bifurcation buckling with the increasing grid point numbers. It is assumed the arch has an included angle of $2\mathrm{\Theta }=90$, common slenderness ratio in practical engineering S/r_x = 100 or 150. As observed, the critical buckling load rapidly approaches stable as the grid point number gradually grows. When the nodal grid points number is greater than 15, any further growth of grid points has no significant impact on the calculation results. Thus $n=15$ is used in subsequent analysis to ensure ideal results can be determined.

Table 1

Table 1. Convergence of HDQ algorithm for solving critical buckling load.

To further validate the accuracy of the presented solutions, finite element analysis is carried out by using commercial software Ansys. In the analysis, the 188-beam element is utilized to model the arch axis, and a total of 220 elements were meshed. The element shape is checked in advance. Control nodes that are fixed along and perpendicular to the axis are tied to the arch feet node with Combine14 element to simulate the linear elastic supports. Riks algorithm is implemented with proper displacement increment to follow the nonlinear buckling path. Three percent of the anti-symmetric buckling displacement which is determined by eigenvalue analysis is introduced to the nonlinear analysis, tracing the bifurcation path. The convergence of analysis is strictly controlled. The critical buckling load, together with the corresponding buckling evolution process are recorded.

Figure 3 illustrates the comparisons of HDQ and FE results for the dimensionless in-plane critical buckling load $qR/N_{TB}$ against included angle $2\mathrm{\Theta }\in \left[0,90\right]$ in Figure 3a, and in Figure 3b for the dimensionless in-plane critical buckling load against the modified slenderness ratio $\lambda \in \left[0,40\right]$. In these figures, the dash line and chain line represent the upper limit point buckling load and bifurcation buckling load, respectively. FE predictions are marked in hollow circles with red, blue, and orange colors for different slenderness ratios S/r_x = {100, 150, 200}. As observed, the HDQ solutions are in good agreement with FE predictions for shallow arches with $2\mathrm{\Theta }<60$, but it slightly overestimates the critical buckling load for deep arches, with a maximum difference of 3.75%. Hence it can concluded that the nonlinear buckling behavior of circular arches with linear elastic supports can be precisely predicted by the presented HDQ-based solution strategy.

Figure 3

Figure 3. Comparisons with FE results for (a) Buckling load against central angle, (b) Buckling load against modified slenderness ratio.

5 Conclusion

This paper investigates the nonlinear in-plane buckling of shallow steel arches with linear elastic supports. A HDQ-based solution strategy is proposed to follow the complete buckling path. Convergence analysis and FE comparison results confirmed that the presented method has good attributes of low computational cost and high accuracy in handling strong nonlinear problems. The buckling behaviors of arches is clarified. It is found that, under uniformly distributed radial load, the circular arch having linear elastic supports could buckle either in a limit point buckling mode or a bifurcation mode, depending on the modified slenderness ratio. The location of bifurcation point on the primary equilibrium path may occur before or after the upper limit point. Since the shear deformation effect is not considered in the analysis, the current calculation method may overestimate the critical buckling load of circular arches with small slendnerness ratio and large included angle. Future endeavor will be made on the buckling behaviors analysis by employing a new strain form that considers shear deformation to improve the calculation accuracy.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

ZZ: Conceptualization, Methodology, Project administration, Supervision, Writing – original draft, Writing – review and editing. WW: Formal Analysis, Investigation, Software, Validation, Writing – original draft. YL: Conceptualization, Formal Analysis, Funding acquisition, Methodology, Resources, Validation, Visualization, Writing – review and editing. JP: Data curation, Formal Analysis, Funding acquisition, Validation, Writing – original draft. WX: Data curation, Funding acquisition, Methodology, Resources, Writing – review and editing.

Funding

The authors declare that financial support was received for the research and/or publication of this article. This paper is financially supported by National Natural Science Foundation of China (No. 52508168 & No. 52279127), Guangdong Basic and Applied Basic Research Foundation (No. 2022A1515111135 & No. 2025A1515010407), and Youth Innovation Talent Project of Ordinary Colleges and Universities in Guangdong Province (No. 2025KQNCX233).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Wang CM, Zhang JM. Optimal shape of an arch under a central point load for maximum buckling load. J Eng Mech (2024) 150:06024002. doi:10.1061/JENMDT.EMENG-7750

CrossRef Full Text | Google Scholar

2. Machacek J. Buckling lengths of steel circular arches respecting non-uniform arch axial forces. Thin-Walled Structures (2022) 180:109916. doi:10.1016/j.tws.2022.109916

CrossRef Full Text | Google Scholar

3. Pi Y-L, Bradford MA. Non-linear buckling and postbuckling analysis of arches with unequal rotational end restraints under a central concentrated load. Int J Sol Structures (2012) 49:3762–73. doi:10.1016/j.ijsolstr.2012.08.012

CrossRef Full Text | Google Scholar

4. Pi Y-L, Bradford MA. Nonlinear analysis and buckling of shallow arches with unequal rotational end restraints. Eng Structures (2013) 46:615–30. doi:10.1016/j.engstruct.2012.08.008

CrossRef Full Text | Google Scholar

5. Han S-C, Ham H-D, Kanok-Nukulchai W. Geometrically non-linear analysis of arbitrary elastic supported plates and shells using an element-based Lagrangian shell element. Int J Non-linear Mech (2008) 43:53–64. doi:10.1016/j.ijnonlinmec.2007.09.011

CrossRef Full Text | Google Scholar

6. Plaut RH. Buckling of shallow arches with supports that stiffen when compressed. J Eng Mech (1990) 116:973–6. doi:10.1061/(ASCE)0733-9399(1990)116:4(973)

CrossRef Full Text | Google Scholar

7. Cai J, Feng J. Buckling of parabolic shallow arches when support stiffens under compression. Mech Res Commun (2010) 37:467–71. doi:10.1016/j.mechrescom.2010.05.004

CrossRef Full Text | Google Scholar

8. Zhang Q, Cai J, Xia W, Feng J. Snap-through of shallow circular arches with variable horizontal supports under unilateral displacement control. J Eng Mech (2020) 146:04020013. doi:10.1061/(ASCE)EM.1943-7889.0001747

CrossRef Full Text | Google Scholar

9. Zhou Y, Yi Z, Stanciulescu I. Nonlinear buckling and postbuckling of shallow arches with vertical elastic supports. J Appl Mech (2019) 86:061001. doi:10.1115/1.4042572

CrossRef Full Text | Google Scholar

10. Lu Y, Cheng Y, Han Q. Experimental investigation into the in-plane buckling and ultimate resistance of circular steel arches with elastic horizontal and rotational end restraints. Thin-Walled Structures (2017) 118:164–80. doi:10.1016/j.tws.2017.05.010

CrossRef Full Text | Google Scholar

11. Han Q, Cheng Y, Lu Y, Li T, Lu P. Nonlinear buckling analysis of shallow arches with elastic horizontal supports. Thin-Walled Structures (2016) 109:88–102. doi:10.1016/j.tws.2016.09.016

CrossRef Full Text | Google Scholar

12. Hu C-F, Li Z, Hu Q-S. On non-linear behavior and buckling of arch-beam structures. Eng Structures (2021) 239:112214. doi:10.1016/j.engstruct.2021.112214

CrossRef Full Text | Google Scholar

13. Hu C-F, Huang Y-M. In-plane nonlinear elastic stability of pin-ended parabolic multi-span continuous arches. Eng Structures (2019) 190:435–46. doi:10.1016/j.engstruct.2019.04.013

CrossRef Full Text | Google Scholar

14. Zhang Z, Liu A, Fu J, Yang J, Liu Y, Huang Y. A theoretical study on nonlinear in-plane buckling of shallow angle-ply laminated arches with elastic supports. Compos Structures (2021) 269:114009. doi:10.1016/j.compstruct.2021.114009

CrossRef Full Text | Google Scholar

15. Kiss LP, Messaoudi A. Assessments of the non-linear instability of arches with imperfect geometry. Structures (2025) 71:108031. doi:10.1016/j.istruc.2024.108031

CrossRef Full Text | Google Scholar

16. Zhang Z, Liu Y, Liu L, Liu A, Chen Z, Yang X. Dynamic buckling of functionally graded porous graphene platelet reinforced composite arches under a locally distributed radial load. Thin-Walled Structures (2025) 208:112838. doi:10.1016/j.tws.2024.112838

CrossRef Full Text | Google Scholar

17. Shagholani Loor A, Rabani Bidgoli M, Hamid M. Optimization and buckling of rupture building beams reinforced by steel fibers on the basis of adaptive improved harmony search-harmonic differential quadrature methods. Case Stud Construction Mater (2021) 15:e00647. doi:10.1016/j.cscm.2021.e00647

CrossRef Full Text | Google Scholar

18. Mosallaie Barzoki AA, Ghorbanpour Arani A, Kolahchi R, Mozdianfard MR, Loghman A. Nonlinear buckling response of embedded piezoelectric cylindrical shell reinforced with BNNT under electro–thermo-mechanical loadings using HDQM. Composites B: Eng (2013) 44:722–7. doi:10.1016/j.compositesb.2012.01.052

CrossRef Full Text | Google Scholar

19. Safarpour M, Safarpour H, Civalek O. Wave propagation analysis of functionally graded bio-composite circular plates using an improved sinusoidal shear deformation theory resting on an advanced viscoelastic foundation. Eur J Mech - A/Solids (2025) 112:105688. doi:10.1016/j.euromechsol.2025.105688

CrossRef Full Text | Google Scholar

20. Tong G, Pi Y-L, Bradford MA, Tin-Loi F. In-Plane nonlinear buckling analysis of deep circular arches incorporating transverse stresses. J Eng Mech (2008) 134:362–73. doi:10.1061/(ASCE)0733-9399(2008)134:5(362)

CrossRef Full Text | Google Scholar