Abstract

This study aimed at obtaining a semianalytical solution for nonlinear dynamic system of shallow arches. Taylor method was applied to find the analytical solution, and an investigation of their dynamic characteristic was carried out to verify the applicability of this methodology for the shallow arches under step or periodic excitation. A polynomial solution can be obtained from this multistep approach with respect to time, and direct buckling as well as indirect buckling of the shallow arches can be observed, also. The results indicated that the dynamic buckling load level was higher with higher shape factor. Additionally, a change of attractor in phase space was investigated. Coupling in symmetric mode as well as asymmetric mode was observed in case of indirect buckling, and a sensitive response was also manifested during sinusoidal and beating excitation. These results of applying multistep Taylor series for the investigation of displacement response and attractor change revealed that this analytical approach was valid in explaining the dynamic buckling behavior of shallow arches under direct and indirect snapping.

1. Introduction

Large span roof structures are generally designed as arches and spherical shells. In other words, the structural performance of the arch depends on its shape. Thus, lighter, thinner, and larger space can be made with less materials compared to traditional flat structures. Shallow arches have received considerable attention in the field of architecture because of their beautiful shape and economic efficiency. However, the response to transverse loading of the shallow arch is quite different from that of the flat beam. The structural design also requires dynamic behavior analysis under various loads. Since a large roof can be treated separately from the substructures, the arch’s boundary is usually designed by hinge, and with shallow arch roofs the primary interest is in the instability of vertical direction behavior. The dynamic unstable behavior of shallow arches is generally complex depending on the initial condition, e.g., geometrically imperfection, and shape of the arches. In particular, direct snapping in symmetric mode and indirect snapping due to coupling in asymmetric mode manifest very sensitive behavior and progress to chaotic behavior, attracting deep interest from many researchers.

Investigation of dynamic instability of shallow arches has started with the fundamental study of Hoff and Bruce [1] in the mid-1950s and was followed by many researches afterwards in the beginning (Humphreys [2]; Lock [3]; Hsu [4]; Ariaratnam and Sankar [5]; Sundararajan and Kumani [6]; Donaldson and Plaut [7]; and Kounadis et al. [8]). These early researches on shallow arches resulted in rebounding studies of chaotic motion or global dynamic behavior (Blair et al. [9]; Levitas et al. [10]), accurate solution of free-vibration (Tseng et al. [11]; De Rosa and Franciosi [12]), internal resonance (Bi and Dai [13]; Lacarbonara and Rega [14]), and buckling under boundary condition and moving load (Kong et al. [15]; Chen and Lin [16, 17]). Recently, many investigations for design and construction of shallow arch roofs have been conducted, including the identification and stability for shallow arches (Ha et al. [18]), and sensitivity of shallow arch (Virgin et al. [19]). On the other hand, the initially curved or buckled microbeams in mechanics (Farokhi and Ghayesh [20]) have been studied in terms of nonlinearity, geometrical imperfection, resonance, and chaotic behaviors. Various numerical researches have been conducted on the nonlinear behaviors of the microbeams and their application. However, this buckled configuration of the microbeams is due to some external effects such as temperature increase or compressive forces, which are not associated with the arch roof. Even though in the study of arch roofs the range of arch shapes and sizes differs from that of microbeams, these investigations are still useful for all engineering fields. The evaluation of dynamic buckling of shallow arch roof generally uses energy method or phase space of system approach or direct approach (Lin and Chen [21]) of investigating significant growth of displacement response (Budiansky and Roth [22]) by obtaining dynamic response of the governing equations in a wide range of parameters. This process of dealing with dynamic loading requires numerical method (Belytschko [23]; Argyris et al. [24]) or analytical techniques for solving nonlinear differential equations, and the process of formulating the equations differs by the methodology. One of the direct purposes of dynamic analysis is finding out the transient response to external force and ultimately attaining the solution to motion equations. However, it is difficult to obtain an analytical solution for dynamic system of shallow arches under the various external forces.

The techniques to obtain an analytical solution to nonlinear differential equations have been considerably investigated in the past (Sadighi et al. [25]). Analytical or semianalytical methods include the traditional Taylor’s power series method (Barrio [26]; Barrio et al. [27]), Adomian decomposition method (Adomian and Rach [28]; Adomian and Rach [29]), homotopy analysis method (Liao [30]), and homotopy perturbation method (He [31]; He [32]; and Chowdhury et al. [33]). Unlike numerical ones, these methods are limited in their applicability and are heavily dependent on the influence of their parameters. The Taylor series method, which has the longest history, provides an accurate analytical solution with excellent stability. Recently, Barrio et al. [27] explained that the Taylor series method provides advantages including easy formation with the method of using variable order and step-size and usefulness in many dynamic systems requiring a high-precision solution. In particular, when we need to perform long-time numerical simulations and sometimes to look for high precision solutions of differential systems, the family of ordinary differential equation solvers that can answer the requirements is the Taylor Series Method (Barrio et al. [34]). Additionally, it has been reported that appropriately adjusting the number of terms and error limit is advantageous in computing the error limit and attaining an accurate solution fairly easily (Barrio et al. [27]; Shon et al. [35]). Traditional Taylor series method requires an infinite series expansion for attaining periodic response of a dynamic system and can be inappropriate for discrete nonlinear governing equations. However, the Taylor series method expanding in multistep can compute a high precision analytical solution with relatively less degree of differential coefficients. Therefore, this study applied multistep Taylor series method to compute an accurate semianalytical solution of a dynamic governing equation of motion for shallow sinusoidal arches and examined their dynamic response and attractor of phase space to determine its applicability. This paper is organized as follows. Sections 2 and 3 describe the formulation of nonlinear governing equation of motion for shallow arches and theoretical analysis technique of multistep Taylor series method, respectively. Section 4 computes a multistep Taylor series solution for the nonlinear governing equation of motion under symmetric mode and investigates the dynamic behavior of shallow arches. Section 5 examines the dynamic response and attractor in phase space of the system for an indirect snapping model under asymmetric mode and periodic excitation. Last, Section 6 concludes this paper.

2. Nonlinear Equations of Motion

The governing equation of motion is expressed in the following equation when the shape and coordinate system of the shallow arches are defined as in Figure 1 (Lin and Chen [21]; Ha et al. [18]). In this figure, the two ends of the arch are separated from each other by the distance . The original shape of the unloaded arch is denoted by , and the vertical and longitudinal deflections are expressed as a function of and , respectively. The arch system is subjected to vertical loading within the coordinate system, and a dimensionless parameter, , is used. Besides, the constant parameters of the system are density of mass , Young’s modulus , cross-sectional area , and moment of inertia .

Figure 1 

The shape of a pin-ended shallow sinusoidal arch under distributed dynamic loading.

The kinetic energy of the shallow arches can be expressed as follows:The potential energy (strain energy) of the system consists of two parts, as shown in (2), one () due to the axial force and one () due to bending. The energy can be expressed as follows:where axial strain and curvature of the shallow arches using Green’s strain tensor are as follows:There are two types of external forces exerting on the arch: excitation force and the damping force. Dissipative force, damping, is considered a uniformly distributed viscous damping force only in the transversal direction. As such, the external works by the external force and dissipative force are as follows:where is the viscous damping coefficient. Now, consider that axial force is constant , and and the boundary condition at , is . Then, the can be expressed as follows:The variations of potential energy , kinetic energy , and work can be formulated as follows:Equations (7)-(9) are substituted to extended Hamilton’s principle given bywhich yields the following nonlinear partial differential equations governing the motion of the shallow arch system (Ha et al. [18]):The following parameters are used to express the above equation (11) with dimensionless variables. Here, represents the radius of gyration of the cross-section.Thus, (11) above can be expressed as the following by using dimensionless system of Figure 1.The initial shape, , deflection , and transverse loading of sinusoidal arches are defined as the following. The representing assumed deformed shape of the arch satisfies the boundary condition (; ) of a pin-ended arch.Accordingly, substitution of the assumed shape, deflection, and load to equation (13) results in nonlinear governing equations of motion for the coefficients of the dimensionless deflection function.where is Kronecker delta, and solving the above equation can explain the dynamic behavior of the arch. Additionally, since and are amplitudes of symmetric mode and asymmetric mode, respectively, the indirect snapping phenomenon can be explained by considering the term.

3. Taylor Series Using Order and Step Parameter

This chapter discusses the formulation of Taylor series analytical method and its governing equations. Equation (15) is defined so as to explain the Taylor series method for computing the solution to initial value problem of nonlinear governing equations and its approximate Taylor series solution. represents an open interval including and is a sufficiently smooth function on . Then, is an analytic function that can be expressed by the following equation.The remainder term for the analytical series solution up to n-degree term is explained and expressed as the following equation, and , which satisfies the equation below, exists in the open interval. Moreover, the error is defined as expressed in the following equation (21).The Taylor series method defines the solution to (15) as n-degree series of (19) except for term, and each differential coefficient, , can be obtained from (15) to solve the equation. The computed solution has the precision with the error range of (21) in the defined interval of . It is very efficient to solve the equation in multistep with and of (21) as the parameters. Adjusting for and is referred to variable order (VO) scheme and variable step-size scheme, respectively. In particular, the two parameters and can determine the error limit, and the accuracy and computational speed of the solution are determined by them. This study defines the solution to (15) as finite Taylor series expressed in (22), and multistep approach is used for computing the solution.The initial value of step in the above equation can be found from previous step, and the differential coefficients, , can be obtained from (15). Accordingly, an approximate analytical solution can be computed from the Taylor series at each step, and the number of terms increases with higher DOF (Degrees of Freedom) and higher order. This study used sufficient order and step-size to compute an accurate solution, and the process of (21) is not further discussed in this paper.

4. Dynamic Snapping Model under Symmetric Mode

Generally, the initiation of snapping process of shallow arches is differentiated by two different mechanisms. That is, one is direct snapping under symmetric mode as illustrated in Figure 2(a), and the other is indirect snapping under asymmetric mode as illustrated in Figure 2(b). In the figure, and are dimensionless initial configuration and dimensionless deflection from the initial configuration, respectively. First of all, symmetric deformation under sinusoidal distributed load is considered to deal with the solution for the case of direct snapping, and the coupling under asymmetric mode is discussed in the next chapter. Accordingly, when the 1st mode of deflection function is considered, (15) is developed into 2nd-order differential equation with respect to as expressed in the following equation.When the dimensionless time parameter is denoted by and the initial values are and , the first two terms of the right-hand side become and . Moreover, when the initial values are substituted into (23), the remaining differential coefficients up to n-degree term can be computed. The following equations show the coefficients of 3rd and 4th terms.Substituting these coefficients to (22) gives us the solution to the polynomial with respect to time, and the number of terms increases with degrees of freedom and order.

(a)

(b)

(a)
(b)

Figure 2 

Deflection modes of the shallow arches: (a) symmetric mode (1st mode) and (b) asymmetric mode (2nd mode).

The appropriateness of the series solution to (23), which is composed of differential coefficients for explaining the dynamic snapping of shallow arches, is evaluated by applying step excitation. The solution was obtained by applying order n and step-size of 7 and 0.1, respectively, to the multistep Taylor series, and the time duration of the computation was set to be 100.

Firstly, the analysis result with (undamped system) is examined. Figure 3(a) shows the analysis result of step excitation with shape parameter of 5, and it shows the displacement response of the shallow arches with the load level increasing from 6 to 16 in increment of 2. Figures 3(a) and 3(b) are the results of the multistep Taylor series method and of the 4th-order Runge-Kutta method (RK4), respectively, for comparing the analysis results. In this case, the time parameters of the Taylor series method and RK4 were set at the same . As shown in the figure, the results of the Taylor series method agree with the results of RK4. Besides, the figure shows that the nonlinear response is clearly differentiated in two groups of load levels (one group of , 8.0, and 10.0 and the other group of , 14.0, and 16.0). The first group had the increase in amplitude of the displacement and gradual lengthening of the period with the increase in the load level, but the second group exhibited conspicuous increase in displacement in comparison with the first group as well as shortening of the period. Moreover, the amplitude of the displacement did not vary linearly with the increase in , and this is the effect of geometrical nonlinearity of the shallow arches. In other words, the former group represents prebuckling load level, and the latter group composes postbuckling load level.

(a)

(b)

(a)
(b)

Figure 3 

Displacement response of the shallow arches subjected to step excitation: (a) multistep Taylor series method and (b) 4th-order Runge-Kutta method (RK4).

Secondly, the maximum displacement response with , 3, 5, and 7 was examined in order to evaluate the solution to the dynamic behavior of the specimen in response to various shape parameter values, and the result is shown in Figure 4. In this figure, dynamic buckling load levels of , 5, and 7 correspond to , 11.05, and 27.9, respectively. In particular, the figure shows that dynamic buckling (snapping) does not take place for the case of , and the boundary point for the occurrence of buckling is higher with higher .

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 4 

Maximum displacement curves under step excitation: (a) ; (b) ; (c) ; and (d) .

The attractor in phase plan as shown in Figure 5 was such that the orbit changed from the form of one drop of water to the form of two drops of water linked together after reaching dynamic snapping load level, and this well explains the dynamic unstable phenomenon of the shallow arches.

(a)

(b)

(a)
(b)

Figure 5 

Phase diagram of the vertical motion for shape parameter : (a) trajectory of prebuckling load level, ; and (b) trajectory of postbuckling load level, .

To explain why snapping does not take place for the case of , there is a need to investigate the equilibrium points of the shallow arch system, and their stability. Let and . Equation (23) can then be written as the 1st-order system, as shown in (26). For this system, the equilibrium points are obtained by solving the system and , and the stability of the points can be observed based on the characteristic polynomial of the Jacobian matrix.In the case of this system, the equilibrium points and their stability were studied by Ha et al. [18], and they can be summarized as follows.

The problem of finding the equilibrium points can be reduced to the investigation of the roots of the 3rd-order polynomial equation as follows:Let , and there are two cases depending on . The first case is for . The equation has only one real root, and the sign is equal to . The second case is for . Then for . The equation has a local maximum and minimum , as shown in (28), and the number of equilibria and their signs depend on . For instance, it has a positive root and two negative roots for .From analyzing the characteristic polynomial of the Jacobian matrix, the stability of the equilibrium point , where is a root of (26), can be defined as follows:(i)Case : unique equilibrium is asymptotically stable.(ii)Case : for , equilibrium is asymptotically stable. For , there are two asymptotically stable points satisfying , and one unstable point satisfying . For , there are two equilibrium points: stable point and asymptotically stable point with .

As regards the condition of and , which means acting on the negatively vertical direction, Figure 6 shows the stability of the equilibrium points under this condition. If and , then the equilibrium point is the unique root of and it is asymptotically stable. In spite of the increase in the magnitude of , the point is still asymptotically stable (see Figure 6(a)). If and , the stability depends on the relationship between and . There are three cases: , , and . For the first case, the equilibrium point under excitation started an asymptotically stable point (see Figure 6(b)). For the next case, a stable point and an asymptotically stable point are observed (see Figure 6(c)). For the final case, two asymptotically stable points and an unstable point are started (see Figure 6(d)). As the increased for these cases, the stabilities are varied, as shown in Figures 6(b)–6(d). As shown in Figure 6, the case of has only a unique asymptotically stable point, and it can be expected that snapping will not take place. On the other hand, the case of started two asymptotically stable points and one unstable point, and it can be expected that the system is very sensitive.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 6 

Equilibrium points under the step excitations: (a) ; (b) ; (c) ; and (d) .

Next, the result of considering the effect of damping coefficient is shown in Figure 7, and it shows that the displacement converges faster with greater . The attraction shown in Figure 7(b) deviates from the trajectory of the first form of a drop of water and converges gradually to the water drop at the opposite side. This indicates that the trajectory becomes stable state asymptotically.  In this model, the shape parameter and load level are 5 and 11.3, respectively, with damping coefficients of , 0.02, and 0.05.

(a)

(b)

(a)
(b)

Figure 7 

Results of the dynamic analysis in consideration of damping coefficient under step excitation : (a) transient response of the motion and (b) trajectory in phase plan.

As the arches may violently vibrate around critical buckling level and at low , there is a need to investigate the equilibrium values about , which is plotted in Figure 8. Figures 8(a) and 8(b) are for and , respectively, and the equilibrium values are computed. The boundaries of this map are for and for , respectively. In the figures, the red points correspond to , and the blue points correspond to . The figures show that the sign of is clearly separated along the critical values when > 0.1, and there appears a mixed critical boundary around for and around for . These results show that the lower damping coefficient is, the more sensitive the variation around the critical boundary will be.

(a)

(b)

(a)
(b)

Figure 8 

Equilibrium points and critical buckling boundary in accordance with the change in the damping coefficient: (a) and (b) .

Finally, consider the stability of the arches in accordance with the different initial values. In particular, the equilibrium values about , which is plotted in Figure 9, are investigated. is the initial value of the 1st-order system, as shown in (26). The damping oscillation with and is analyzed, and the equilibrium values with , 3, 5, and 7 are computed. In the figure, plane element is divided into elements as the initial values on the lattice point of the coordinate and time .

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 9 

Attraction of the 1st-order system of arches under excitation level and damping coefficient : (a) ; (b) ; (c) ; and (d) .

As shown in Figure 6, the equilibrium points are distinguished by the range of . For and , there is a unique stable point, and when the time increases from the initial value , focuses on a single point, as shown in Figures 9(a) and 9(b). For and , there are two asymptotically stable points and an unstable point, and when the time increases from the initial value , focuses on the two stable points, as shown in Figures 9(c) and 9(d). If the trajectory is attracted to the equilibrium points, the lattice element of the initial value is drawn in red, otherwise, in blue. For Figures 9(a) and 9(b), red is almost equal to blue, but Figures 9(c) and 9(d) have different signs and magnitudes.

For better understanding, the extended phase diagrams are drawn in Figure 10. These results can be easily observed from Figures 10(c) and 10(d). From the results of and , the boundary of the color becomes the boundary of the attractor. In particular, a pair of insets can be clearly seen through the equilibrium values, and these points are divided into positive and negative values.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 10 

Expended phase space of the 1st-order system of arches: (a) ; (b) ; (c) ; and (d) .

In conclusion, the analytical series solution converged even when damping was considered, and it could well predict the phenomenon and characteristics of the expected direct snapping of shallow arches.

5. Dynamic Snapping Model under Asymmetric Mode

This chapter describes indirect snapping model due to coupling in symmetric mode and asymmetric mode. Since the influence of asymmetric mode as shown in Figure 2(b) should consider both the first term and second term of displacement function, the following 2nd-order differential equations with respect to and are derived from (15). When and the initial values are set to be , , , and as described in the previous chapter, the first two terms of the power series solution are and for the case of and and for the case of . The remaining terms are expressed in the following when the order is 4 as mentioned in the previous chapter.The analysis continued to investigate the result of the solution to indirect snapping for the shape of with parameters of and . In this analysis, 0.1% of h was applied to account for the initial imperfection, and beating excitation and step and periodic excitation was carried out.

Figure 11 is the displacement response to step excitation for load level and shape parameter . In case of Figure 11(a), imperfection is not considered, but it is introduced by , i.e., 0.1% of , in case of Figure 11(b). Figure 11(b) in consideration of initial imperfection shows that displacement rapidly changes due to the influence of coupling under asymmetric mode past the time level of 10. This is also observed in Figure 12, showing the change of attraction in the phase plan. Considering the boundary of dynamic critical load level for direct snapping, , the result of Figure 11 indicates that the shallow arch test specimens are very sensitive to initial condition.

(a)

(b)

(a)
(b)

Figure 11 

Transient response of shallow arches with initial imperfection () (MTSM): (a) perfection, ; and (b) initial imperfection .

(a)

(b)

(a)
(b)

Figure 12 

Phase diagram of the dynamic response of the analytical model with initial imperfection (, ) (MTSM): (a) trajectory of and (b) trajectory of .

To verify and compare with the analysis results by multistep Taylor series method, the models in Figure 11 are recomputed using RK4, and the results are shown in Figures 13 and 14. In the case of the perfect model, the time parameters of the Taylor series method and RK4 were set at the same . As shown in Figures 11(a) and 13(a), the results of the Taylor series method agree with the results of RK4. Besides, the trajectory in Figure 12(a) well matches the trajectory in Figure 14(a).

(a)

(b)

(c)

(a)
(b)
(c)

Figure 13 

Transient response of the shallow arches with initial imperfection (RK4): (a) perfection , ; (b) initial imperfection , ; and (c) initial imperfection , .

(a)

(b)

(c)

(a)
(b)
(c)

Figure 14 

Phase diagram of the dynamic response of the analytical model with initial imperfection (, ) (RK4): (a) perfection , ; (b) initial imperfection , ; and (c) initial imperfection , .

Considering the imperfection model in Figure 11(b), the results of RK4 corresponding to this imperfection model are shown in Figures 13(b) and 13(c), which are for and for , respectively. From the result shown in Figure 13(b), the response pattern is similar, but it can be seen that they do not coincide. On the other hand, Figure 13(c) is more consistent than Figure 13(b). As the imperfection model is more sensitive than the perfection model, a small time internal is required for RK4. Besides, the results pattern can be easily observed from Figures 14(b) and 14(c). The results shown in Figures 13 and 14 indicate that the shallow arches are very sensitive to the initial condition, and higher accuracy is required to analyze the sensitive model.

Multistep Taylor series were applied to the analysis results of sinusoidal excitation and beating excitation as shown in Figure 15. The periodic parameters and of 1.0 and 0.05, respectively, were applied as shown in Figure 15.In (34) and (35), is the first natural angular frequency of the model, and the excitation is applied at the same period as the natural frequency of the case of .