A Semianalytical Approach for Nonlinear Dynamic System of Shallow Arches Using Higher Order Multistep Taylor Method

1Department of Architectural Engineering, Korea University of Technology and Education, Cheonan 31253, Republic of Korea 2Acrovision Co. Ltd., B-1411, 70 Dusan-ro, Geumcheon-gu, Seoul 08584, Republic of Korea 3Interdisciplinary Program in Creative Engineering, Korea University of Technology and Education, Cheonan 31253, Republic of Korea 4School of Liberal Arts, Korea University of Technology and Education, Cheonan 31252, Republic of Korea


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.

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  0 (), 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 ().
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  0 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 (=  0 ), and / = 0 and the boundary condition at  = 0,  is (0) = () = 0.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 by which 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.
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  1 and  2 are amplitudes of symmetric mode and asymmetric mode, respectively, the indirect snapping phenomenon can be explained by considering the  2 term.

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  0 and  is a sufficiently smooth function on  ×   ×   .Then, () is an analytic function that can be expressed by the following equation.
The remainder term (, ,  0 ) 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).
(, ,  0 )     ≤ 1 ( + 1)! max The Taylor series method defines the solution to (15) as ndegree series of ( 19) except for  term, and each differential coefficient,  () ( 0 ), 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 [ 0 ,  0 +  ℎ ].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 stepsize 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,  () ( 0 ), 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.

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,  0 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  1 as expressed in the following equation.
When the dimensionless time parameter  is denoted by  0 and the initial values are  1 ( 0 ) =  0 and Ḋ1 ( 0 ) = Ḋ 0 , the first two terms of the right-hand side become  0 and Ḋ0 ( −  0 ).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.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  = 0 (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  ℎ = 0.1.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 Λ = 6.0, 8.0, and 10.0 and the other group of Λ = 12.0, 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.
Secondly, the maximum displacement response with ℎ = 1, 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 ℎ = 3, 5, and 7 correspond to Λ = 3.20, 11.05, and 27.9, respectively.In particular, the figure shows that dynamic buckling (snapping) does not take place for the case of ℎ = 1, and the boundary point for the occurrence of buckling is higher with higher ℎ.
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.
To explain why snapping does not take place for the case of ℎ = 1, there is a need to investigate the equilibrium points of the shallow arch system, and their stability.Let  1 =  − ℎ and u = V.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 V = 0 and V = 0, 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 () =    0 = √(ℎ 2 − 4)/3.The equation has a local maximum (− 0 ) and minimum ( 0 ), 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 0 < 4(ℎ + Λ) < (− 0 ).
As regards the condition of ℎ > 0 and Λ ≤ 0, which means acting on the negatively vertical direction, Figure 6 shows the stability of the equilibrium points under this condition.If ℎ ≤ 2 and Λ = 0, then the equilibrium point is the unique root of () = 4ℎ 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 ℎ > 2 and Λ = 0, the stability depends on the relationship between 4ℎ and (− 0 ).There are three cases: 2 < ℎ < 4, ℎ = 4, and ℎ > 4. For the first case, the equilibrium point under excitation Λ = 0 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 ℎ ≤ 2 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 ℎ > 4 started two asymptotically stable points and one unstable point, and it can be expected that the system is very sensitive.
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.0, 0.02, and 0.05.
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 These results show that the lower damping coefficient  is, the more sensitive the variation around the critical boundary will be. Finally, consider the stability of the arches in accordance with the different initial values.In particular, the equilibrium values about ( =0 , V =0 ), which is plotted in Figure 9, are investigated.( =0 , V =0 ) is the initial value of the 1st-order system, as shown in (26).The damping oscillation with  = 0.As shown in Figure 6, the equilibrium points are distinguished by the range of ℎ.For ℎ < 4 and Λ = 0, there is a unique stable point, and when the time increases from the initial value ( =0 , V =0 ), () focuses on a single point, as shown in Figures 9(a) and 9(b).For ℎ > 4 and Λ = 0, there are two asymptotically stable points and an unstable point, and when the time increases from the initial value ( =0 , V =0 ), () 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 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 ℎ = 5 and ℎ = 7, 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.
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.

Dynamic Snapping Model under
Asymmetric Mode Asymptotically stable influence of asymmetric mode as shown in Figure 2(b) should consider both the first term  1 and second term  2 of displacement function, the following 2nd-order differential equations with respect to  1 and  2 are derived from (15).
When  =  0 and the initial values are set to be  1 ( 0 ) =  0 1 , Ḋ1 ( 0 ) = Ḋ0 1 ,  2 ( 0 ) =  0 2 , and Ḋ2 ( 0 ) = Ḋ 0 2 as described in the previous chapter, the first two terms of the power series solution are  0 1 and Ḋ 0 1 ( −  0 ) for the case of  1 and  0 2 and Ḋ 0 2 ( −  0 ) for the case of  2 .The remaining terms are expressed in the following when the order  is 4 as mentioned in the previous chapter. 1 (3) Mathematical Problems in Engineering The analysis continued to investigate the result of the solution to indirect snapping for the shape of ℎ with parameters of  = 7 and  ℎ = 0.1.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 Λ = 18.5 and shape parameter ℎ = 7.In case of Figure 11(a), imperfection is not considered, but it is introduced by  2 = 0.001ℎ, 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, Λ = 27.9, the result of Figure 11 indicates that the shallow arch test specimens are very sensitive to initial condition.
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   and 13(c), which are for  ℎ = 0.1 and for  ℎ = 0.001, 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),  0 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  = 1.0.This study considered the case of  = 1.0 only.Accordingly, since resonance was expected by excitation, the analysis continued with lower level Λ 0 of 8.0 than the buckling load level under step excitation.Figure 16 shows the analysis result.In this example, the shape parameter and load level of sinusoidal excitation are ℎ = 7 and Λ 0 = 8.0, respectively.Imperfection is not considered in Figure 16(a), but it is introduced by  2 = 0.001ℎ, i.e., 0.1% of ℎ in case of Figure 16(b).After the time elapses to 20 as shown in figure (b), the amplitude of symmetric mode  1 is increased by the perturbation of the asymmetric mode  2 .It indicates that, for the case of sinusoidal excitation without the consideration of imperfection (Figure 16(a)), there was no rapid increase in the observed displacement.Nevertheless, the effect of coupling was very great for the case of considering asymmetric mode (Figure 16(b)).This result is manifested clearer with the trajectory in phase plan as shown in Figure 17.In other words, although there was no change of trajectory in the phase plan as shown in Figure 17(a), Figure 17(b) shows that the  1 curve changes the form of trajectory after the elapse of a certain time.
The analysis of beating excitation was carried out with ℎ = 7 like the case of sinusoidal excitation, and Λ 0 of 8.0, the same load level as for sinusoidal excitation, was applied.Figure 18 shows the analysis result.After the time elapses Analyses of dynamic response to step, sinusoidal, and beating excitation under indirect snapping condition and the influence of asymmetric mode were carried out in order to evaluate the effectiveness of the multistep Taylor series method.Coupling under nonlinearity and symmetric and asymmetric modes was well manifested in dynamic response to these excitations.The dynamic response under periodic excitation at the same period as the natural frequency  0 was very reliable and verified.

Conclusion
This paper examined a dynamic analysis of nonlinear governing equations for such structures as shallow arches and the computational method for their analytical solution.The Taylor method was used to obtain a semianalytical solution to a motion equation, and a solution composed of a polynomial function with respect to time was formulated in multistep.A dynamic analysis was conducted to evaluate the computed solution for the behavior of shallow arches subjected to step and periodic excitations.Besides, to verify the results using the multistep Taylor series method, the 4th-order Runge-Kutta method was used to compute the responses and trajectories of the shallow arches under the excitation, and the results are well matched.The analysis result for the arch under symmetric mode revealed that the dynamic buckling load level became higher with the higher shape parameter value of the arch, and the phase diagram showed a change of attraction at the load level of direct snapping compared to the level prior to the occurrence of direct snapping.The stability of the symmetric mode arches is also summarized, and the equilibria according to ℎ are investigated.For ℎ ≤ 2, dynamic snapping does not happen, and the arches system for ℎ > 4 are very sensitive.In particular, for ℎ > 4, it can be clearly seen that the attractors are divided into positive and negative values.For the case of indirect buckling, the similar attractor and response, which manifested the coupling phenomenon under symmetric and asymmetric modes, were also obtained.A response result sensitive to the initial imperfection was observed in the tests of the sinusoidal and beating excitations,

Figure 1 :
Figure 1: The shape of a pin-ended shallow sinusoidal arch under distributed dynamic loading.
) and 9(b), red is almost equal to blue, but Figures 9(c) and 9(d) have different signs and magnitudes.

Figure 7 :Figure 8 :
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.

Figures 13 and 14 .
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  ℎ = 0.1.As shown in Figures11(a) and 13(a), the results of the Taylor series method agree with the results of RK4.Besides, the trajectory in Figure12(a) well matches the trajectory in Figure14(a).Considering the imperfection model in Figure11(b), the results of RK4 corresponding to this imperfection model are shown in Figures13(b) and 13(c), which are for  ℎ = 0.1 and for  ℎ = 0.001, respectively.From the result shown in Figure13(b), the response pattern is similar, but it can be seen that they do not coincide.On the other hand, Figure13(c) is more consistent than Figure13(b).As the imperfection model 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  ℎ = 0.1.As shown in Figures11(a) and 13(a), the results of the Taylor series method agree with the results of RK4.Besides, the trajectory in Figure12(a) well matches the trajectory in Figure14(a).Considering the imperfection model in Figure11(b), the results of RK4 corresponding to this imperfection model are shown in Figures13(b) and 13(c), which are for  ℎ = 0.1 and for  ℎ = 0.001, respectively.From the result shown in Figure13(b), the response pattern is similar, but it can be seen that they do not coincide.On the other hand, Figure13(c) is more consistent than Figure13(b).As the imperfection model