A Novel Shape Finding Method for the Main Cable of Suspension Bridge Using Nonlinear Finite Element Approach

The determination of the final cable shape under the self-weight of the suspension bridge enables its safe construction and operation. Most existing studies solve the cable shape segment-by-segment in the Lagrangian coordinate system. This paper develops a novel shape finding method for the main cable of suspension bridge using nonlinear finite element approach with Eulerian description. The governing differential equations for a three-dimensional spatial main cable is developed before a one-dimensional linear shape function is introduced to solve the cable shape utilizing the Newton iteration method. The proposed method can be readily reduced to solve the two-dimensional parallel cable shape. Two iteration layers are required for the proposed method. The shape finding process has no need for the information of the cable material or cross section using the present technique. The commonly used segmental catenary method is compared with the present method using three cases study, i.e., a 1666-m-main-span earth-anchored suspension bridge with 2D parallel and 3D spatial main cables as well as a 300-m-main-span self-anchored suspension bridge with 3D spatial main cables. Numerical studies and iteration results show that the proposed shape finding technique is sufficiently accurate and operationally convenient to achieve the target configuration of the main cable.


Introduction
The advancements of the fundamental theory, construction technique, and new materials have enabled the rapid growth of the bridge main span length to cross the large rivers, wide canyons, and deep straits [1]. As a bridge type with largest span capacity, the suspension bridge is always a preferable candidate for a single span longer than one kilometer from the perspectives of the mechanical property and economic performance. Though the longest suspension bridge is still the Akashi Kaikyō Bridge with main span of 1991 m built in 1998, several super-long suspension bridges with the main span longer than two kilometers are recently being constructed, such as the 2023-m-main-span Çanakkale 1915 Bridge in Turkey, the 2018-m-main-span Shiziyang Bridge in China, and the 2300-m-main-span Zhanggao Bridge across Yangtze River in China. The gravity stiffness of the main cable owing to the self-weight or initial tension stress is one of the most important reasons for the continuous growth of the main span for the suspension bridge. Its proportion in the stiffness of the suspension bridge system will also increase with the growth of the main span length.
Essentially, the gravity stiffness steams from the geometric nonlinearity, which strongly depends on the configuration and internal force of the main cable. However, the shape of the main cable should be prescribed accounting for the effects of all external forces, which is always unable to be amended during the construction [2]. The configuration or target shape for the main cable at the initial equilibrium state under dead loads due to the self-weights of the bridge should be predetermined at the preliminary design stage. A shape finding or form finding process is therefore conducted to determine the shape and internal forces of the main cable through minimizing the dead-load-induced deformation of the bridge. Moreover, the increase of the main cable dimension and the reduction of the cable design safety factor require a more refined estimation of the cable configuration and internal forces in the structural components [3].
In the past several decades, the theory of the shape finding for the main cable of suspension bridge has received intensive attention. For a given cable segment of specified unstrained length, cross section area, Young's elastic modulus, and density, there are many kinds of relations between node force and node displacement/position, e.g., a link element considering sag effect [4,5], a two-node curved cable element [6,7], and various two-dimensional [8] and three-dimensional catenary cable elements [9][10][11][12]. The segmental parabola theory assumes the shape of the main cable to be a parabolic curve under uniform load [4], [13]. This approximation ignores the concentrated forces transferred from the hangers, which is only acceptable for medium and small span suspension bridges.
The segmental catenary theory allows one to achieve the analytical solutions of each cable segment shape. The equations describing the analytical relations between the axial forces and the strained/unstrained lengths of a cable segment under the action of the selfweight can be found in many pioneering studies, e.g., [14]. These equations were then widely employed to develop various shape finding approaches for the main cable of the suspension bridge, such as the initial force method or segmental catenary method (SCM) [3,[15][16][17][18], the target configuration under dead load (TCUD) method [19], the improved TCUD method [20], the Generalized TCUD method [21], the coordinate iteration method [22], and the perturbation approach [23,24]. The unstrained length of each main cable segment in these methods is unknown and solved in the successive nonlinear equations using the nonlinear finite element iteration [15,16,[25][26][27][28][29]. Generally, at least two loop layers are required for calculating the unstrained length and axial force of each catenary segment before determining the shape of the main cable. To avoid a loop over each element, the unstrained length of each cable is treated as an unknown parameter in formulating a tangential stiffness matrix in the TCUD method [19] and the Generalized TCUD method [21]. Shape finding process based on the analytical solution of each catenary element is mostly exact. However, its theoretical derivation is complex, and the achievement of the algorithm requires high programing skills. Moreover, most of them are restricted to two-dimensional (2D) parallel main cables.
The three-dimensional (3D) spatial main cables, which are intended to enhance the lateral stiffness of the suspension bridge, have been employed in some recent projects, such as the Yongjong Grand Bridge in Korea, Jiangdong Bridge in China, the New Oakland Bay Bridge in the USA, and Halogaland Bridge in Norway. Tang et al. [25] extended the 2D segmental catenary approach to 3D by iterating the axial component of the main cable to achieve a target sag-to-span ratio. Kim et al. [20] estimated the spatial cable shape of the Yongjong Grand Bridge using an improved TCUD method coupled with the initial force method by introducing a 3D elastic catenary cable element. The procedures in these studies are lengthier and more time-consuming for spatially curved cables as compared with the planar cable system.
To simplify the calculation procedure and reduce the computation burden, the finite element method (FEM) was employed to estimate the 3D main cable shape. A five-step algorithm was proposed in Xiao et al. [30] and Song et al. [31] using the common FEM software ANSYS to solve both 2D planar and 3D spatial cable shapes, which is convenient to be executed at the preliminary design stage. However, the use of the common FEM program results in the insufficiently theoretical analysis for the core algorithm.
As can be seen, most above-mentioned methods are developed in the Lagrangian coordinate system, which expresses the quantities of interest as the function of the undeformed position or unstrained length of the main cable segment. However, the unstrained length of each cable segment is usually unknown, resulting in more loop layers. In this paper, a shape finding procedure for the main cable using the NFEM and the Eulerian description is developed. A modified governing equation for the whole main cable, rather than cable segments, based on the force equilibrium is proposed by introducing a condition function to express the concentrated forces from hangers [14]. All quantities of the main cable are expressed in a fixed coordinate system. That is, if the longitudinal direction along the bridge axis, the vertical and lateral directions of a suspension bridge are defined as the x-, y-and z-axes, the variables of the main cable related to the y and z coordinates of the cable are expressed by a fixed x coordinate instead of the unstrained length of the cable. An FEM discretization process [32] starting from a linear basic function is then conducted to solve the governing partial differential equations to get a group of nonlinear algebraic equations. The iteration schemes for solving the cable shape with known and unknown horizontal axial force are proposed. The commonly used segmental catenary method is compared with the present approach. The shapes of 2D parallel and 3D spatial main cables, which are applied to a 1666-m main span earth-anchored suspension bridge are solved using the present approach and the segmental catenary method, respectively. The self-anchored Yongjong Grand Bridge analyzed by Kim et al. [16] is also examined using the proposed method to validate its accuracy. The proposed algorithm provides an alternative to estimate the main cable shape at the preliminary design stage using only two-layer iteration and less programming. The use of the FEM allows the present method to be readily embedded in some commonly used FE analysis software and easily used by general bridge engineers.

Governing Equations of Spatial Main Cable
As shown in Figure 1a, a spatial main cable is described in a Eulerian coordinate system. The x-axis is the longitudinal direction along the bridge axis. The y-axis is in the vertical direction and the z-axis is in the lateral direction of the bridge. The tension forces of the cable at two ends are decomposed into three components, i.e., , , and , , . The unit weight of the main cable is defined as q. The projected length of the cable in the x-direction is L. The total number of hangers is n. The tension of ith (i = 1, 2, …, n) inclined hanger is noted as , which consists of a vertical force of in the ydirection and a lateral component of in the z direction, as shown in Figure 1b At the final state of the bridge, the lumped force transferred from each inclined hanger can be determined from the dead load of the girder and self-weight of the hanger [25]. The distribution force due to the gravity of the main cable are also prescribed. To establish the equilibrium equations of a spatial main cable, four assumptions are employed in this study: (1) The main cable of the suspension bridge is ideally flexible and unable to bear any bending moment. Only the axial tension is considered for the main cable. (2) The axial deformation of the main cable is tiny such that the area of the cross section for the main cable remains unchanged and the cable weight per unit length is a constant. (3) The linear elastic assumption for all materials is adopted. (4) All the loadings on the main cable are parallel to the y-o-z plane. The tension in the x-direction of the main cable is constant and defined as .
The first three assumptions are commonly employed during the shape finding process at the preliminary design stage, e.g., Li et al. [12]; Kim et al. [16,20]; Zhang et al. [18]. As reported by Zhang et al. [18], the flexural stiffness of the main cable is insignificant and can be neglected during the construction process. Moreover, the main cable should work substantially below the yield stress such that the axial deformation is tiny and the variation of the cross-section area is negligible. For most modern suspension bridges, the hangers in are designed to incline only in the cross-sectional plane of the bridge. There is no or insignificant tension component of each hanger along the longitudinal direction. Therefore, the tension of the main cable in the x-direction remains constant.
The tension force along the axis of the main cable is assumed as , which can be decomposed into the , and in the x-, y-and z-direction, respectively. The direction of is parallel to the cable element at the position of interest, i.e., = , , . The ratio between the , and can be expressed as: Taking the derivative of Equations (1) and (2) with the respect to gives the form of The distributions of and along the x-direction can also be formulated as the functions of real loadings: in which − is a condition function defined as: It is noteworthy that the Einstein summation convention is employed hereafter in this study for the purpose of the formula simplification, e.g., − = ∑ − in Equations (5) and (6). The governing equations for the spatial main cable can be obtained by combing the Equations (3) and (6) with the form of in which the dash denotes spatial differentiation with respect to . For the two-dimensional main cable that is parallel to the x-o-y plane, a constant z coordinate is assigned.

Nonlinear FEM Solution for the Cable Shape
A differential equation with the strong form is defined as , , , , ⋯ = 0, in which ∎ is a function symbol, is the independent variable, is the dependent variable, , , ⋯ are the first-, second-and high-order derivatives with respect to the . In FEM, it can be discretized using a basis function or a shape function to produce a weak form of ∫ , , , , ⋯ Ω = 0, in which Ω is the computational domain of interest. To solve the cable shape, a one-dimensional linear shape function is adopted in this study with the form of = + (12) in which and are undetermined coefficients estimated by the nodal coordinates of the element. By using the Galerkin's method, the cable shape in terms of and can be formulated as: By introducing the , the discrete algebraic form or weak form of Equations (9) and (10) can be expressed as: in which the subscript (j = 1, 2, …, n) indicates the hanger ID, which is introduced to distinguish with the subscript before the Einstein summation convention is applied.
The subscripts k and l in the following sections are also utilized for the same purpose. The integration by parts is applied to Equations (15) and (16), which leads to in which BC indicates the integral boundary condition. The substitution of Equation (8) into Equations (17) and (18) produces: in which and are the final y-direction and z-direction coordinates of the main cable. An additional formula for describing the ratio between and is given as: As formulated in Equation (12), the differential of a one-dimensional linear shape function with respect to x is a constant for an element. Meanwhile, the node of the main cable that is located at the top of the pylon is assumed to be fixed. The Dirichlet boundary condition is therefore applied to the variables of y and z. The integral form of Equations (19) and (20) can be alternatively solved using a form of the summation as: in which dL is the length of element Ω in the x-direction. It is worth mentioning that ( ) and ( ) are constant within each element. However, their values are not identical for different elements. To solve this nonlinear equation related to and , the Newton iteration method [33] is introduced: in which ( ) is the function to be solved, is the approximate solution of ( ) = 0 at (n + 1)-step iteration, is a relaxation coefficient, which is set as 1.0 in this study. ′( ) in Equations (22) and (23) is a matrix for a multivariable problem, which consists of the partial derivatives of the objective function with respect to each unknow variables: in which and are functions described by the first members of Equations (22) and (23), which are derived from the equilibrium equations in y-direction (Equation (9) The nodal coordinates of the main cable are readily updated using the summation of the ( ) and ( ). It is worth noting that he initial nodal coordinates ( , ) before the iteration can be assigned as arbitrary values, the coordinates of the intersection point (IP) between the pylon and main cable are suggested for the operational convenience. The above iteration scheme is illustrated in Figure 2 with a given , which is called the inner-layer iteration in this study.

However, in most cases,
is an undetermined parameter. The middle span-sag or the coordinate of the middle point in y-direction of the main cable is always prescribed. An outer-layer iteration is therefore implemented, as described in Figure 2. At each outer-layer iteration, the coordinate of the middle point in y-direction can be solved by a given using the inner-layer iteration. A one-to-one mapping relationship = ( ) is assumed, in which (∎) is the function symbol. As a result, the solution of can be determined by finding the root of the equation ( ) = ( ) − = 0. Since ( ) has no explicit expression, the Secant method is utilized to achieve a numerical solution of . The initial conditions for and are suggested as the 1.0 time and 1.1 times of For two-dimensional parallel main cables, the coordinate of a main cable in z-direction is an invariant constant. The can reduced to a × matrix as = , in which is assembled from the by removing the -related terms in Equation (25) with the form of Figure 2. Flowchart of the iteration scheme for solving the cable shape.

Comparison with the Segmental Catenary Method
The segmental catenary method is an exact solution for an element of the main cable deformed due to its self-weight. The lengths of each segment in x-y-and z-directions are solved as: in which , and are the distances between the two nodes of a segment along the global x-y-and z-axes, , and are three components of the global nodal force, is the elastic modulus of the main cable, is the cross-sectional area, is the self-weight of the main cable per unit length, is the unstrained length of the cable segment. Figure 3 illustrates the iteration scheme of the segmental catenary method. For a two-pylon suspension bridge, a guess of the axial force = , , at one end of the main cable is made first. The unstrained length and and of first segment are obtained by solving the transcendental equations described in Equations (34) and (36). The axial force of next segment, which will be used to solve the and and is then readily calculated by introducing the tensile force of the first hanger.
will be adjusted by comparing the coordinates of the IP at the other pylon in y-and z-directions, i.e., and with the calculated and ̂ as well as the with to achieve a convergent solution. A comparison between the SCM and present NFEM is shown in Table 1. The SCM focuses on each cable segment to develop its analytical expression of the cable shape using algebraic equations (Equations (34) and (36)) in the Lagrangian coordinate system with respect the unstrained cable length . The present NFEM solves the differential equations of the whole cable in the Eulerian coordinate system with respect the bridge axis coordinate x using the nonlinear finite element approximation. The SCM requires the inputs of the axis stiffness of the cable, i.e., to determine the unstrained cable length while no information regarding the cable material or cross section is needed in the NFEM for solving the cable shape. It is noteworthy that the unstrained cable length, which is of great importance during the construction stage, can also be calculated using the NFEM by subtracting the elongation from the final length of the cable.

Case 1: 2D Parallel Main Cables of an Earth-Anchored Suspension Bridge
An earth-anchored suspension bridge with 2D parallel main cables is utilized, as shown in Figure 4.   Two-layer iterations are implemented to perform the shape finding process for the main cable by following the scheme in Figure 2. The outer-layer and inner-layer iteration tolerances are set as = 1.0 × 10 and = 1.0 × 10 . Figure 6 illustrates the iteration results, in which a convergent horizontal component of the axial force, i.e., is achieved in five outer-layer iteration steps. Between each adjacent outer-layer iterations, less than four-step inner-layer iterations are required for solving the cable shape. The main cable shape from the initial input coordinates of to the final shape finding result that corresponds the final outer-and inner-layer iteration is plotted in Figure 7. The shape finding using the commonly used segment catenary method (SCM) is also conducted as a benchmark for comparison purpose. Both the final main cable shapes obtained from this study and the SCM are given in Figure 8a, showing a reasonable agreement with each other. The maximum relative difference at all hanger nodes is 0.0011%. Moreover, the parabolic curve of the main cable in the middle span calculated by the = is also compared with the present solution, as shown in Figure 8b. The difference between these two methods is also insignificant with the maximum relative error of 0.0948%. This is mainly because of the sag-span ratio of this suspension bridge about equal to 1/10 such that the dead loadings along the along the bridge approximately follow the uniform distribution. Furthermore, the static displacement of the final main cable under the action of self-weight and hanger tension forces is analyzed using the FEM, as shown in Figure 9. As can be seen, all nodal displacements of the main cable are less than 0.007 m. It is sufficiently small to believe that the main cable has achieved the target configuration. It is noteworthy that the time costs in this case for both the NFEM and the SCM are less than one second.

Case 2: 3D Spatial Main Cables of an Earth-Anchored Suspension Bridge
Three dimensional spatial main cables are applied to the same suspension bridge used in Case 1. The spacings between two cables at the position of the anchorage and the IP points in z-direction are set as 63.0 m and 3.0 m, respectively, as shown in Figure 10. The vertical (y-direction) force components of all hangers given in Figure 5 are identical to Case 1. Same tolerances are also employed in this case, say = 1.0 × 10 and = 1.0 × 10 . The iteration process is shown in Figure 11. Similar to the Case 1, five outerlayer iterations are implemented to obtain a convergent . All inner-layer iterations are less than five steps to meet the tolerance requirement. The main cable shape from the initial input coordinates in y-and z-directions to the final target configurations is illustrated in Figure 12. The final shape of the 2D parallel main cables in Case 1 is also plotted in Figure 12a. As can be seen, the shapes in y-direction of 2D planar cable and 3D spatial cable match well with each other. This can be explained by comparing the governing equations for the cable shape described by Equations (9) and (11). The difference of the cable shapes in x-o-y plane for parallel and spatial cases in mainly attributed to the . For the long-span suspension bridge studied in this case, the slope in z-direction of the spatial cable is small as compared to the slope in y-direction, resulting in that the effect of the in Equations (9) and (10) is insignificant. This can also be demonstrated from the final main shape in Figure 12. The calculated maximum and using the central difference method are 0.2593 and 0.0033, respectively.   Figure 12. The main cable shape from the initial coordinates to the target configuration for 3D spatial main cables: (a) y-direction; (b) z-direction. Figure 13 illustrates the solution of hanger tension components in z-direction along the bridge span. The positive tension indicates that the inclination of the hanger tilts inward. The maximum positive tension in z-direction occurs at the middle span with the value of 331.637 kN. The negative tensions for outward inclined hangers at two side spans reach the maximum for first or end hanger. This is because the vertical component of the tension for first or end hanger as shown in Figure 5 is large. Figure 14 compares the cable shape between this study and the results obtained by SCM. The maximum relative differences at all hanger nodes in y-and z-direction are 0.0022% and 0.7318%, respectively, suggesting a reasonable agreement with each other. To further validate the accuracy of the final cable shape, a finite element model of the spatial main cable is also established. The static displacement under the action of self-weight and hanger tension forces is illustrated in Figure 15. It can be found that the maximum displacements in y-and z-direction of the main cable are 0.008 m and 0.012 m, which are sufficiently small and acceptable for the engineering applications.  Another case is the Yongjong Grand Bridge in Korea, whose shape finding was investigated by many pioneering studies [16,22,26]. The Yongjong Grand Bridge is the first highway and railway self-anchored suspension bridge with three dimensional main cables. The total length of the bridge is 550 m with a 300-m main span and two 125-m side spans. The sag in middle span is 60 m.

Case 3: 3D Spatial Main Cables of a Self-Anchored Suspension Bridge
has a uniform value of 3.048 KN. The crosssectional area of the main cable is 0.1355 m 2 . The y and z coordinate values of the IP point are = 114.573 m and = 1.50 m, respectively by setting the origin of the Eulerian coordinate system at the midpoint of the bridge with zero elevation. The main cable is divided into 24 segments by 23 hangers. The coordinates of the connection points between hangers and the beam, i.e., and are estimated from the data given by Kim et al. [16] and Luo et al. [26], as listed in Table 2. The IP point is located at x = 125.0 m, and the mid-point of the main span is at x = 275 m.
The iteration process for the Yongjong Grand Bridge using present method is illustrated in Figure 16. Both the outer-layer and inner-layer iterations meet the tolerance requirement in five steps. The solved horizontal component of the cable axial force is 48,149.1 kN, which is close to the value, 49,541 kN given by Kim et al. [16] with the relative error of 2.89 %. The solution of the cable shape is compared with the results obtained by Kim et al. [16] and Luo et al. [26] in Table 3. As can be seen, the relative differences in ydirection between the present study and Kim et al. [16] are less than 0.027%. In z-direction, the absolute values of the relative difference are less than 0.340%. The location of the IP was further studied by Luo et al. [26], which is different from the coordinates used in the present study and Kim et al. [16]. The maximum absolute value of the relative difference between present study and Luo et al. [26] in y and z directions are 0.041% and 0.524%, respectively, which are slightly larger than the errors between the present study and Kim et al. [16]. However, the present method is proved to be accurate enough to find a reasonable cable shape of the suspension bridge for engineering applications.