A Novel Analytical Model for Structural Analysis of Long-Span Hybrid Cable-Stayed Suspension Bridges

Abstract

The hybrid cable-stayed suspension bridge is used to combine the advantages of cable-stayed and suspension bridges and hence has a broad prospect for application. The conventional simplified analytical models of the hybrid bridge are usually developed based on a schematic with the cable-stayed and suspension systems working separately without any overlapping zone, which cannot represent the modern hybrid bridge system. In this study, a novel analytical model is proposed based on the modified suspension–elastic foundation beam theory to estimate the mechanical performance and deflection of the hybrid bridge system with the consideration of the overlapped section between the suspension and stayed cables. The governing equations of the hybrid bridge system are developed based on the elastic foundation beam theory and the deflection theory, which are derived separately in the hybrid section, the pure suspension section and the cable-stayed section. The general solution of each section is presented. The Transfer Matrix Method is then employed to solve the unknowns from one end to the other, which are in turn used to solve the internal forces of the hybrid bridge system caused by the concentrated load. In addition, in view of no variation in the unstressed length of the main cable, the compatibility equation of the main cable is established with consideration of the longitudinal displacement of the main tower, which is used to derive the formulas for the internal force and deflection of the hybrid system. The model can be easily complied in any programming platform, such as Matlab, with simple input parameters, which can eliminate the complex finite element modeling process. Hence, it can be easily used in the preliminary design stage to determine the optimal size and layout of the bridge. Then, a case study is presented for the verification of the proposed model under a vertical load, which is simplified from the Xihoumen Bridge, a combined highway and railway bridge with a main span of 1488 m. Good agreement is obtained between the proposed model and the finite element method. Meanwhile, it is found that there exists a negative deflection zone for the main beam at a distance from the concentrated vertical load, which is mainly caused by the deflection of the main cables, leading to the cambering of the beam.

1. Introduction

A hybrid bridge with both cable-stayed and suspension systems is a collaboration of a cable-stayed bridge and a suspension bridge, which, with the optimal design, is able to combine the optimal properties of these two bridge systems [1,2,3,4]. The suspension cable system can guarantee stable and safe erection processes, and meanwhile, the additional stayed cables can increase the stiffness and reduce the beam deformability. Hence, the hybrid system is able to improve the span capacity and structural stability and reduce the pylon height and anchorage volume [5].

The hybrid system has a long history of development for about 170 years, with the first system proposed in 1855 by American bridge designer Roebling for the Niagara Falls highway railway bridge [6], which has a main span of 251 m. The first hybrid system is also called the Roebling system where the force is mainly taken by suspenders, while the stay cables are used to increase the stiffness and bearing capacity of the suspension bridge. In 1883, the second hybrid bridge system was proposed by Franz Dichinger, which was named as the Dichinger system [7]. In this system, the cable-stayed bridge part is independent of the suspension bridge part, where the cable-stayed bridge sustains the loads on both sides of the side span and the middle-span, while the suspension bridge takes the load of the middle part of the middle span with a ground-anchored suspension system. With continuous research and development, novel cable-stayed suspension hybrid bridge systems have been proposed, such as the Troitsky system [8], the Lin Tongyan system [9], and the Gimsing system [10]. Meanwhile, the span capacity of the hybrid bridge system has been improved significantly by the design optimization to take full advantage of both cable-stayed and suspension bridges. Currently, the largest-span hybrid bridge is the Third Bosphorus Bridge with a span length of 1408 m [11]. The Xihoumen Rail-cum-Road Bridge, which is currently under construction, has a span of 1488 m [12].

While the hybrid bridge system has a long history of development, only a few cases of applications have been constructed. In addition, its longest span is still far less than that of a suspension bridge, which is the 1915 Canakkale Bridge with a span of 2023 m [13]. The reasons could be the complex force transformation mechanism between various load-carrying components, including the main cables, stayed cables, beam, and tower, where the forces shared by various components are unclear. Meanwhile, a lack of systematic experimental and theoretical analyses, as well as the optimization scheme, could also be another important concern.

For the complex hybrid bridge system, the most straightforward method to determine the mechanical static and dynamic behavior is through the finite element method (FEM), which is used to combine the modeling techniques of cable-stayed and suspension bridges. With proper modeling techniques, the deflection and internal forces of the hybrid system can be solved directly [14,15,16]. However, for the complex hybrid system with a long span, the FEM involves complex modeling elaboration with a significant amount of input data and complex compiling process, which are all susceptible to computational errors. Only the engineers with a sophisticated modeling background and experience can complete the modeling process. In addition, in the preliminary design stage, a significant amount of design work will focus on the optimization of the geometric configuration of the hybrid system, resulting in varying configurations with different geometric parameters, which is not easily adjusted in the FEM.

In order to facilitate the preliminary optimization process of the hybrid system, researchers also make many efforts on the development of an efficient analytical model with simple input parameters based on elastic beam theory and the potential energy method. For example, Hegab [17] proposed an iterative program to estimate the static response of cable-stayed bridges based on the potential energy method. Aboul-ella [18] extended Hegab’s method to consider the flexibility of the supporting towers. Xia and Zhang [19] developed an analytical formula for the gravity stiffness and cable stiffness of the hybrid bridge system based on the energy relation principle associated with deformation of the cable. Thai and Choi [20] proposed an analytical method to estimate the ultimate strength of suspension bridges, which could consider material and geometric nonlinearities. Cheng and Li [21] proposed a simplified method to predict the deflections of the hybrid bridge system under live loads based on the potential energy method. The deflection of the main beam was described using a Fourier series, and a few iteration procedures were required. However, their approach cannot be used for any longitudinal configuration not satisfying the restriction of low displacements in the tower. Wang et al. [22] proposed an iterative form-finding method by integrating the FEM and analytical formulas with optimization algorithms for the target configuration of a new type of spatial self-anchored hybrid cable-stayed suspension bridge under the dead load. Feng et al. [23] put forward an analytical model for the early design stage of cable-stayed suspension bridges based on the Hellinger–Reissner variational method. Zhang et al. [24,25] introduced an analytical algorithm to predict the response of a hybrid-cable-stayed suspension bridge under a transverse live load, based on the assumption that the state of the full bridge under the dead load is known. It is noted that the governing equations in the proposed algorithm are derived and solved based on the conservation of the unstressed length of each section of cable, elevation difference in various spans, and stress balance of the beam. Other related analytical solutions for the hybrid bridge system are referred to in the literature [26,27,28,29,30]. For all aforementioned analytical models, researchers mainly focus on the global deflection of the hybrid system. In addition, for the simplicity of the calculation, the schematic of the hybrid system does not consider the overlapped section between the suspension and cable-stayed parts. In fact, in order to improve the fatigue performance of the end hangers, the modern hybrid bridge system has a long overlapping section.

In this paper, a novel analytical model based on the modified suspension–elastic foundation beam theory is proposed to analyze the mechanical performance and deflection of the hybrid bridge system by considering the overlapping section between the suspension and stayed cables. The governing equations of the hybrid bridge system are derived separately in the hybrid section, the pure suspension section, and the cable-stayed section. Details of the model development and solution derivation are given. A case study is then presented to verify the proposed model using the finite element model. Finally, the characteristics of the influence lines for the vertical deflection, beam rotation, moment, and shear are investigated using the proposed model.

2. Assumptions

Figure 1 shows the elevation of the hybrid bridge system, including the cable-stayed section ($L_{xl}$), the pure suspension section ($L_{xd}$), and the hybrid cable-stayed suspension section ($L_{jc}$). The following assumptions are made in the proposed model.

(a)

The relationships of stress–strain of all materials satisfy Hook’s laws.

(b)

The main cable is in the parabolic curve under the uniform vertical load.

(c)

The hangers and stayed cables are densely distributed, and hence can be treated as a homogeneous membrane.

(d)

The hangers have negligible elongation and inclination under the dead and service loads. Then, the deck deflection is equal to that of the suspended cable.

(e)

The horizontal displacement of the suspended cable is neglected.

3. Differential Equilibrium Equations

Due to different loading systems in different sections of the hybrid bridge, the differential equilibrium equations are derived in the hybrid, suspension, and cable-stayed sections, respectively.

3.1. Hybrid Section

Figure 2 shows the force diagram of the infinitesimal element for the hybrid section. In the figure, N, Q, and M are the axial force, shear force, and moment of the main beam, respectively; $w$ and $\nu$ are the deflection of the main cable and main beam, respectively; q and p are the distributed load caused by the dead and the live load, respectively; c is the distributed force of the stay cable based on the homogeneous membrane assumption of stayed cables; $c_y$ is the vertical component of c; $c_{yd}$ and $c_{yp}$ are the vertical components of the stay cable generated by the dead and live loads, respectively, and $c_y=c_{yd}+c_{yp}$. The horizontal force H of the main cable is composed of two parts: $H_q$ and $H_p$, which are generated by the dead and live loads, respectively, $H=H_q+H_p$.

Taking the bending moment at the left end point of the main beam, it gives the following:

$$H\left(dy+dw\right)+dM-Ndv-\left(V+Q\right)dx=0$$

(1)

where V is the vertical component of the suspension cable. Based on the equilibrium in the vertical direction, it yields the following equation:

$$d\left(V+Q\right)+\left(p+q-c_y\right)dx=0$$

(2)

Meanwhile, for the beam, the differential equation for the deformed shape is as follows:

$$M=-EI{\nu }^{″}$$

(3)

where E is the elastic modulus of the main beam, and I is the bending moment of inertia of the main beam. Taking the first derivative of Equation (1) and substituting Equations (2) and (3) into Equation (1), it gives the following:

$$EI{\nu }^{\left(4\right)}-(H_q+H_p)(y^{″}+w^{″})+N{\nu }^{″}-(p+q-c_y)=0$$

(4)

Based on the assumption that the hangers have negligible elongation, $w=\nu$, and $H_q{y}^{″}=c_{yd}-q$ [31], then, the Equation (4) can be simplified as follows:

$$EI{\nu }^{\left(4\right)}-\left(H_q-N\right){\nu }^{″}-{\mathrm{H}}_{\mathrm{p}}{\nu }^{″}+c_{yp}=p+H_p{y}^{″}$$

(5)

Equation (5) is the differential equilibrium equation of the hybrid section of the bridge system under the assumption of deflection theory. Considering that the main beam of the bridge is supported by the cable in a large area and the area of the pure suspension section in the middle span is generally small, the ratio of the stiffness of the main beam part to the gravity stiffness of the main cable is larger than that of the conventional suspension bridge. In addition, for the bridge of more than 1000 m, $H_p\le 0.2H_q$. Therefore, Equation (5) is further simplified according to the linear deflection theory and the term ${\mathrm{H}}_{\mathrm{p}}{\nu }^{″}$ is ignored.

The vertical component force generated by the stay cable due to the deflection of the main beam is a function of the deflection and position coordinates. To facilitate the solution, it is simplified according to Winkler-type elastic foundation. Assuming the vertical elastic support stiffness provided by the stay cable is k, then, the vertical component force of the stay cable is $c_{yp}=k\nu$. Substituting this into Equation (5) gives the following:

$$EI{\nu }^{\left(4\right)}-\left(H_q-N\right){\nu }^{″}+k\nu =p+H_p{y}^{″}$$

(6)

The general solution of the corresponding linear homogeneous differential Equation (6) with constant coefficients is as follows:

$$\nu =\left\{\begin{array}{c}Ash\left(\gamma x\right)\cos \left(\delta x\right)+Bsh\left(\gamma x\right)\sin \left(\delta x\right)+Cch\left(\gamma x\right)\cos \left(\delta x\right)+Dch\left(\gamma x\right)\sin \left(\delta x\right) H_q<2\sqrt{EIk}\\ Axsh\left(\eta x\right)+Bxch\left(\eta x\right)+Csh\left(\eta x\right)+Dch\left(\eta x\right) H_q=2\sqrt{EIk}\\ Ash\left(\eta x\right)+Bch\left(\eta x\right)+Csh\left(\epsilon x\right)+Dch\left(\epsilon x\right) H_q>2\sqrt{EIk}\end{array}\right.$$

(7)

where

$$\left\{\begin{array}{c}\begin{array}{c}\gamma =\sqrt{{\displaystyle \frac{2\sqrt{EIk}+H_q-N}{4EI}}}, \delta =\sqrt{{\displaystyle \frac{2\sqrt{EIk}-\left(H_q-N\right)}{4EI}}} H_q-N<2\sqrt{EIk}\\ \\ \end{array}\\ \eta =\sqrt{{\displaystyle \frac{H_q-N+\sqrt{(H_q-N{)}^2-4EIk}}{2EI}}}, \epsilon =\sqrt{{\displaystyle \frac{H_q-N-\sqrt{(H_q-N{)}^2-4EIk}}{2EI}}}, H_q-N\ge 2\sqrt{EIk}\end{array}\right.$$

The vertical displacement $\nu$, the first derivative of displacement ${\nu }^{\prime }$, the bending moment M, and the shear force Q of any point of the main beam in the hybrid section form a characteristic mechanical vector $T=(\nu , {\nu }^{\prime }, M,Q{)}^T$, and its relationship with the coefficient vector $S=(A, B, C, D{)}^T$, is as follows:

$$T=K_{jc}S$$

(8)

where $K_{jc}$ is the eigen mechanical matrix.

(1)

When $H_q-N<2\sqrt{EIk}$

$$K_{jc}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right]$$

and the value of each element is as follows:

$$\begin{array}{c}{a}_{11}=sh\left(\gamma x\right)\mathit{cos}(\delta x), a_{12}=sh(\gamma x\left)\mathit{sin}(\delta x\right), a_{13}=ch\left(\gamma x\right)\mathit{cos}(\delta x), a_{14}=ch(\gamma x\left)\mathit{sin}(\delta x\right)\hfill \\ a_{21}=\gamma ch\left(\gamma x\right)\mathit{cos}(\delta x)-\delta sh(\gamma x\left)\mathit{sin}(\delta x\right), a_{22}=\gamma ch\left(\gamma x\right)\mathit{sin}(\delta x)+\delta sh(\gamma x\left)\mathit{cos}(\delta x\right)\hfill \\ a_{23}=\gamma sh\left(\gamma x\right)\mathit{cos}(\delta x)-\delta ch(\gamma x\left)\mathit{sin}(\delta x\right), a_{24}=\gamma sh\left(\gamma x\right)\mathit{sin}(\delta x)+\delta ch(\gamma x\left)\mathit{cos}(\delta x\right)\hfill \\ a_{31}=EI\left\{\right({\delta }^2-{\gamma }^2\left)sh\right(\gamma x\left)\mathit{cos}(\delta x\right)+2\gamma \delta ch\left(\gamma x\right)\mathit{sin}(\delta x\left)\right\}\hfill \\ a_{32}=EI\left\{\right({\delta }^2-{\gamma }^2\left)sh\right(\gamma x\left)\mathit{sin}(\delta x\right)-2\gamma \delta ch\left(\gamma x\right)\mathit{cos}(\delta x\left)\right\}\hfill \\ a_{33}=EI\left\{\right({\delta }^2-{\gamma }^2\left)ch\right(\gamma x\left)\mathit{cos}(\delta x\right)+2\gamma \delta sh\left(\gamma x\right)\mathit{sin}(\delta x\left)\right\}\hfill \\ a_{34}=EI\left\{\right({\delta }^2-{\gamma }^2\left)ch\right(\gamma x\left)\mathit{sin}(\delta x\right)-2\gamma \delta sh\left(\gamma x\right)\mathit{cos}(\delta x\left)\right\}\hfill \\ a_{41}=EI\left\{\right(3{\delta }^2-{\gamma }^2\left)\gamma ch\right(\gamma x\left)\mathit{cos}(\delta x\right)+(3{\gamma }^2-{\delta }^2)\delta sh\left(\gamma x\right)\mathit{sin}(\delta x\left)\right\}\hfill \\ a_{42}=EI\left\{\right(3{\delta }^2-{\gamma }^2\left)\gamma ch\right(\gamma x\left)\mathit{sin}(\delta x\right)-(3{\gamma }^2-{\delta }^2)\delta sh\left(\gamma x\right)\mathit{cos}(\delta x\left)\right\}\hfill \\ a_{43}=EI\left\{\right(3{\delta }^2-{\gamma }^2\left)\gamma sh\right(\gamma x\left)\mathit{cos}(\delta x\right)+(3{\gamma }^2-{\delta }^2)\delta ch\left(\gamma x\right)\mathit{sin}(\delta x\left)\right\}\hfill \\ a_{44}=EI\left\{\right(3{\delta }^2-{\gamma }^2\left)\gamma sh\right(\gamma x\left)\mathit{sin}(\delta x\right)-(3{\gamma }^2-{\delta }^2)\delta ch\left(\gamma x\right)\mathit{cos}(\delta x\left)\right\}\hfill \end{array}$$

The inverse matrix is $K_{jc}^{-1}=\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right]$, and the value of each element is as follows:

$$\begin{array}{c}{b}_{11}=-{\displaystyle \frac{({\delta }^2-{\gamma }^2)ch\left(\gamma x\right)\mathit{sin}(\delta x)+2\gamma \delta sh(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta }}\hfill \\ b_{21}={\displaystyle \frac{({\delta }^2-{\gamma }^2)ch\left(\gamma x\right)\mathit{cos}(\delta x)-2\gamma \delta sh(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta }}\hfill \\ b_{31}={\displaystyle \frac{({\delta }^2-{\gamma }^2)sh\left(\gamma x\right)\mathit{sin}(\delta x)+2\gamma \delta ch(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta }}\hfill \\ b_{41}=-{\displaystyle \frac{({\delta }^2-{\gamma }^2)sh\left(\gamma x\right)\mathit{cos}(\delta x)-2\gamma \delta ch(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta }}\hfill \\ b_{12}={\displaystyle \frac{(3{\delta }^2-{\gamma }^2)\gamma sh\left(\gamma x\right)\mathit{sin}(\delta x)+(3{\gamma }^2-{\delta }^2\left)\delta ch\right(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)}}\hfill \\ b_{22}=-{\displaystyle \frac{(3{\delta }^2-{\gamma }^2)\gamma sh\left(\gamma x\right)\mathit{cos}(\delta x)-(3{\gamma }^2-{\delta }^2\left)\delta ch\right(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)}}\hfill \\ b_{32}=-{\displaystyle \frac{(3{\delta }^2-{\gamma }^2)\gamma ch\left(\gamma x\right)\mathit{sin}(\delta x)+(3{\gamma }^2-{\delta }^2\left)\delta sh\right(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)}}\hfill \\ b_{42}={\displaystyle \frac{(3{\delta }^2-{\gamma }^2)\gamma ch\left(\gamma x\right)\mathit{cos}(\delta x)-(3{\gamma }^2-{\delta }^2\left)\delta sh\right(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)}}\hfill \\ b_{13}={\displaystyle \frac{ch\left(\gamma x\right)\mathit{sin}(\delta x)}{2\gamma \delta EI}}, b_{23}=-{\displaystyle \frac{ch\left(\gamma x\right)\mathit{cos}(\delta x)}{2\gamma \delta EI}}, b_{33}=-{\displaystyle \frac{sh\left(\gamma x\right)\mathit{sin}(\delta x)}{2\gamma \delta EI}}, b_{43}={\displaystyle \frac{sh\left(\gamma x\right)\mathit{cos}(\delta x)}{2\gamma \delta EI}}\hfill \\ b_{14}={\displaystyle \frac{\delta ch\left(\gamma x\right)\mathit{cos}(\delta x)-\gamma sh(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)EI}}\hfill \\ b_{24}={\displaystyle \frac{\delta ch\left(\gamma x\right)\mathit{sin}(\delta x)+\gamma sh(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)EI}}\hfill \\ b_{34}=-{\displaystyle \frac{\delta sh\left(\gamma x\right)\mathit{cos}(\delta x)-\gamma ch(\gamma x\left)\mathit{sin}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)EI}}\hfill \\ b_{44}=-{\displaystyle \frac{\delta sh\left(\gamma x\right)\mathit{sin}(\delta x)+\gamma ch(\gamma x\left)\mathit{cos}(\delta x\right)}{2\gamma \delta ({\delta }^2+{\gamma }^2)EI}}\hfill \end{array}$$

(2)

When $H_q-N=2\sqrt{EIk}$

$$K_{jc}=\left[\begin{array}{cccc}xsh\left(\eta x\right)& xch\left(\eta x\right)& sh\left(\eta x\right)& ch\left(\eta x\right)\\ sh\left(\eta x\right)+\eta xch\left(\eta x\right)& ch\left(\eta x\right)+\eta xsh\left(\eta x\right)& \eta ch\left(\eta x\right)& \eta sh\left(\eta x\right)\\ -EI\left\{2\eta ch\right(\eta x)+{\eta }^{2}xsh(\eta x\left)\right\}& -EI\left\{2\eta sh\right(\eta x)+{\eta }^{2}xch(\eta x\left)\right\}& -EI{\eta }^{2}sh\left(\eta x\right)& -EI{\eta }^{2}ch\left(\eta x\right)\\ -EI{\eta }^2\left\{3sh\right(\eta x)+\eta xch(\eta x\left)\right\}& -EI{\eta }^2\left\{3ch\right(\eta x)+\eta xsh(\eta x\left)\right\}& -EI{\eta }^{3}ch\left(\eta x\right)& -EI{\eta }^{3}sh\left(\eta x\right)\end{array}\right]$$

The inverse matrix is:

$$K_{jc}^{-1}=\left[\begin{array}{cccc}-{\displaystyle \frac{ch\left(\eta x\right)}{2}}& {\displaystyle \frac{sh\left(\eta x\right)}{2}}& -{\displaystyle \frac{ch\left(\eta x\right)}{2\eta EI}}& {\displaystyle \frac{sh\left(\eta x\right)}{2{\eta }^{2}EI}}\\ {\displaystyle \frac{sh\left(\eta x\right)}{2}}& -{\displaystyle \frac{ch\left(\eta x\right)}{2}}& {\displaystyle \frac{sh\left(\eta x\right)}{2\eta EI}}& -{\displaystyle \frac{ch\left(\eta x\right)}{2{\eta }^{2}EI}}\\ -{\displaystyle \frac{2sh\left(\eta x\right)-\eta xch\left(\eta x\right)}{2}}& {\displaystyle \frac{3ch\left(\eta x\right)-\eta xsh\left(\eta x\right)}{2\eta }}& {\displaystyle \frac{xch\left(\eta x\right)}{2\eta EI}}& {\displaystyle \frac{ch\left(\eta x\right)-\eta xsh\left(\eta x\right)}{2{\eta }^{3}EI}}\\ {\displaystyle \frac{2ch\left(\eta x\right)-\eta xsh\left(\eta x\right)}{2}}& -{\displaystyle \frac{3sh\left(\eta x\right)-\eta xch\left(\eta x\right)}{2\eta }}& -{\displaystyle \frac{xsh\left(\eta x\right)}{2\eta EI}}& -{\displaystyle \frac{sh\left(\eta x\right)-\eta xch\left(\eta x\right)}{2{\eta }^{3}EI}}\end{array}\right]$$

(3)

When $H_q-N>2\sqrt{EIk}$

$$K_{jc}=\left[\begin{array}{cccc}sh\left(\eta x\right)& ch\left(\eta x\right)& sh\left(\epsilon x\right)& ch\left(\epsilon x\right)\\ \eta ch\left(\eta x\right)& \eta sh\left(\eta x\right)& \epsilon ch\left(\epsilon x\right)& \epsilon sh\left(\epsilon x\right)\\ -EI{\eta }^{2}sh\left(\eta x\right)& -EI{\eta }^{2}ch\left(\eta x\right)& -EI{\epsilon }^{2}sh\left(\epsilon x\right)& -EI{\epsilon }^{2}ch\left(\epsilon x\right)\\ -EI{\eta }^{3}ch\left(\eta x\right)& -EI{\eta }^{3}sh\left(\eta x\right)& -EI{\epsilon }^{3}ch\left(\epsilon x\right)& -EI{\epsilon }^{3}sh\left(\epsilon x\right)\end{array}\right]$$

The inverse matrix is

$$K_{jc}^{-1}=\left[\begin{array}{cccc}-{\epsilon }^{2}sh\left(\eta x\right)& {\displaystyle \frac{{\epsilon }^{2}ch\left(\eta x\right)}{\eta }}& -{\displaystyle \frac{sh\left(\eta x\right)}{EI}}& {\displaystyle \frac{ch\left(\eta x\right)}{\eta EI}}\\ {\epsilon }^{2}ch\left(\eta x\right)& -{\displaystyle \frac{{\epsilon }^{2}sh\left(\eta x\right)}{\eta }}& {\displaystyle \frac{ch\left(\eta x\right)}{EI}}& -{\displaystyle \frac{sh\left(\eta x\right)}{\eta EI}}\\ {\eta }^{2}sh\left(\epsilon x\right)& -{\displaystyle \frac{{\eta }^{2}ch\left(\epsilon x\right)}{\epsilon }}& {\displaystyle \frac{sh\left(\epsilon x\right)}{EI}}& -{\displaystyle \frac{ch\left(\epsilon x\right)}{\epsilon EI}}\\ -{\eta }^{2}ch\left(\epsilon x\right)& {\displaystyle \frac{{\eta }^{2}sh\left(\epsilon x\right)}{\epsilon }}& -{\displaystyle \frac{ch\left(\epsilon x\right)}{EI}}& {\displaystyle \frac{sh\left(\epsilon x\right)}{\epsilon EI}}\end{array}\right]/({\epsilon }^2-{\eta }^2)$$

3.2. Pure Suspension Section

Figure 3 shows the force diagram of the infinitesimal element for the pure suspension section. When compared with the hybrid section in Figure 2, no stay cable force is presented in the pure suspension section.

Since the axial force of the beam is small in the pure suspension section, it is neglected when deriving the differential equation. In addition, no stay cable force is presented in the pure suspension section. Therefore, the differential equation of the pure suspension section can be simplified from Equation (6), and it gives the following:

$$EI{\nu }^{\left(4\right)}-H_q{\nu }^{″}=p+H_p{y}^{″}$$

(9)

It is seen that Equation (9) is consistent with the differential equation of the suspension bridge.

The general solution of the linear homogeneous differential Equation (9) with constant coefficients is as follows:

$$\nu =Ash\left(\alpha x\right)+Bch\left(\alpha x\right)+Cx+D$$

(10)

where

$$\alpha =\sqrt{H_q/EI}$$

(11)

Then, the relationship of the characteristic mechanical vector T and the coefficient vector S is as follows:

$$T=K_{xd}S$$

(12)

where $K_{xd}=\left[\begin{array}{cccc}sh\left(\alpha x\right)& ch\left(\alpha x\right)& x& 1\\ \alpha ch\left(\alpha x\right)& \alpha sh\left(\alpha x\right)& 1& 0\\ -EI{\alpha }^{2}sh\left(\alpha x\right)& -EI{\alpha }^{2}ch\left(\alpha x\right)& 0& 0\\ -EI{\alpha }^{3}ch\left(\alpha x\right)& -EI{\alpha }^{3}sh\left(\alpha x\right)& 0& 0\end{array}\right]$, and its inverse matrix is as follows:

$$K_{xd}^{-1}=\left[\begin{array}{cccc}0& 0& {\displaystyle \frac{sh\left(\alpha x\right)}{EI{\alpha }^2}}& -{\displaystyle \frac{ch\left(\alpha x\right)}{EI{\alpha }^3}}\\ 0& 0& -{\displaystyle \frac{ch\left(\alpha x\right)}{EI{\alpha }^2}}& {\displaystyle \frac{sh\left(\alpha x\right)}{EI{\alpha }^3}}\\ 0& 1& 0& {\displaystyle \frac{1}{EI{\alpha }^2}}\\ 1& -x& {\displaystyle \frac{1}{EI{\alpha }^2}}& -{\displaystyle \frac{x}{EI{\alpha }^2}}\end{array}\right]$$

3.3. Cable-Stayed Section

Figure 4 shows the force diagram of the infinitesimal element for the cable-stayed section. When compared with the hybrid section in Figure 2, there is no contribution from the main cable in the cable-stayed section.

The derivation of the equilibrium differential equation of the cable-stayed section can be derived from the Equation (4) by removing the main cable-related terms and meanwhile substituting $q-c_{yd}=0$. It gives the following:

$$EI{\nu }^{\left(4\right)}+N{\nu }^{″}-(p-c_{yp})=0$$

(13)

Similarly, assuming that the cable is equivalent to a beam supported on an elastic foundation and the vertical elastic support stiffness k provided by the cable, then, the Equation (13) can be simplified as follows:

$$EI{\nu }^{\left(4\right)}+N{\nu }^{″}+k\nu =p$$

(14)

The general solution of the linear homogeneous differential Equation (14) with constant coefficients is as follows:

$$\nu =\left\{\begin{array}{c}Ash\left(\beta x\right)\cos \left(\lambda x\right)+Bsh\left(\beta x\right)\sin \left(\lambda x\right)+Cch\left(\beta x\right)\cos \left(\lambda x\right)+Dch\left(\beta x\right)\sin \left(\lambda x\right) N<2\sqrt{EIk}\\ Axsin\left(\varphi x\right)+Bxcos\left(\varphi x\right)+Csin\left(\varphi x\right)+Dcos\left(\varphi x\right) N=2\sqrt{EIk}\\ Asin\left(\varphi x\right)+Bcos\left(\varphi x\right)+Csin\left(\psi x\right)+Dcos\left(\psi x\right) N>2\sqrt{EIk}\end{array}\right.$$

(15)

where

$$\left\{\begin{array}{c}\begin{array}{c}\beta =\sqrt{{\displaystyle \frac{2\sqrt{EIk}-N}{4EI}}}, \lambda =\sqrt{{\displaystyle \frac{2\sqrt{EIk}+N}{4EI}}} N<2\sqrt{EIk}\\ \\ \end{array}\\ \eta =\sqrt{{\displaystyle \frac{N+\sqrt{N^2-4EIk}}{2EI}}}, \epsilon =\sqrt{{\displaystyle \frac{N-\sqrt{N^2-4EIk}}{2EI}}} N\ge 2\sqrt{EIk}\end{array}\right.$$

Then, the relationship of the characteristic mechanical vector T and the coefficient vector S is as follows:

$$T=K_{xl}S$$

(16)

where $K_{xl}$ is the eigen mechanical matrix for the cable-stayed section.

(1)

When $N<2\sqrt{EIk}$

$K_{xl}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right]$ and each element is computed as follows:

$$\begin{array}{c}{a}_{11}=sh\left(\beta x\right)\mathit{cos}(\lambda x), a_{12}=sh(\beta x\left)\mathit{sin}(\lambda x\right), a_{13}=ch\left(\beta x\right)\mathit{cos}(\lambda x), a_{14}=ch(\beta x\left)\mathit{sin}(\lambda x\right)\hfill \\ a_{21}=\beta ch\left(\beta x\right)\mathit{cos}(\lambda x)-\lambda sh(\beta x\left)\mathit{sin}(\lambda x\right), a_{22}=\beta ch\left(\beta x\right)\mathit{sin}(\lambda x)+\lambda sh(\beta x\left)\mathit{cos}(\lambda x\right)\hfill \\ a_{23}=\beta sh\left(\beta x\right)\mathit{cos}(\lambda x)-\lambda ch(\beta x\left)\mathit{sin}(\lambda x\right), a_{24}=\beta sh\left(\beta x\right)\mathit{sin}(\lambda x)+\lambda ch(\beta x\left)\mathit{cos}(\lambda x\right)\hfill \\ a_{31}=EI\left\{\right({\lambda }^2-{\beta }^2\left)sh\right(\beta x\left)\mathit{cos}(\lambda x\right)+2\beta \lambda ch\left(\beta x\right)\mathit{sin}(\lambda x\left)\right\}\hfill \\ a_{32}=EI\left\{\right({\lambda }^2-{\beta }^2\left)sh\right(\beta x\left)\mathit{sin}(\lambda x\right)-2\beta \lambda ch\left(\beta x\right)\mathit{cos}(\lambda x\left)\right\}\hfill \\ a_{33}=EI\left\{\right({\lambda }^2-{\beta }^2\left)ch\right(\beta x\left)\mathit{cos}(\lambda x\right)+2\beta \lambda sh\left(\beta x\right)\mathit{sin}(\lambda x\left)\right\}\hfill \\ a_{34}=EI\left\{\right({\lambda }^2-{\beta }^2\left)ch\right(\beta x\left)\mathit{sin}(\lambda x\right)-2\beta \lambda sh\left(\beta x\right)\mathit{cos}(\lambda x\left)\right\}\hfill \\ a_{41}=EI\left\{\right(3{\lambda }^2-{\beta }^2\left)\beta ch\right(\beta x\left)\mathit{cos}(\lambda x\right)+(3{\beta }^2-{\lambda }^2)\lambda sh\left(\beta x\right)\mathit{sin}(\lambda x\left)\right\}\hfill \\ a_{42}=EI\left\{\right(3{\lambda }^2-{\beta }^2\left)\beta ch\right(\beta x\left)\mathit{sin}(\lambda x\right)-(3{\beta }^2-{\lambda }^2)\lambda sh\left(\beta x\right)\mathit{cos}(\lambda x\left)\right\}\hfill \\ a_{43}=EI\left\{\right(3{\lambda }^2-{\beta }^2\left)\beta sh\right(\beta x\left)\mathit{cos}(\lambda x\right)+(3{\beta }^2-{\lambda }^2)\lambda ch\left(\beta x\right)\mathit{sin}(\lambda x\left)\right\}\hfill \\ a_{44}=EI\left\{\right(3{\lambda }^2-{\beta }^2\left)\beta sh\right(\beta x\left)\mathit{sin}(\lambda x\right)-(3{\beta }^2-{\lambda }^2)\lambda ch\left(\beta x\right)\mathit{cos}(\lambda x\left)\right\}\hfill \end{array}$$

Its inverse matrix is $K_{xl}^{-1}=\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right]$, and the value of each element is as follows:

$$\begin{array}{c}{b}_{11}=-{\displaystyle \frac{({\lambda }^2-{\beta }^2)ch\left(\beta x\right)\mathit{sin}(\lambda x)+2\beta \lambda sh(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda }}\hfill \\ b_{21}={\displaystyle \frac{({\lambda }^2-{\beta }^2)ch\left(\beta x\right)\mathit{cos}(\lambda x)-2\beta \lambda sh(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda }}\hfill \\ b_{31}={\displaystyle \frac{({\lambda }^2-{\beta }^2)sh\left(\beta x\right)\mathit{sin}(\lambda x)+2\beta \lambda ch(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda }}\hfill \\ b_{41}=-{\displaystyle \frac{({\lambda }^2-{\beta }^2)sh\left(\beta x\right)\mathit{cos}(\lambda x)-2\beta \lambda ch(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda }}\hfill \\ b_{12}={\displaystyle \frac{(3{\lambda }^2-{\beta }^2)\beta sh\left(\beta x\right)\mathit{sin}(\lambda x)+(3{\beta }^2-{\lambda }^2\left)\lambda ch\right(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)}}\hfill \\ b_{22}=-{\displaystyle \frac{(3{\lambda }^2-{\beta }^2)\beta sh\left(\beta x\right)\mathit{cos}(\lambda x)-(3{\beta }^2-{\lambda }^2\left)\lambda ch\right(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)}}\hfill \\ b_{32}=-{\displaystyle \frac{(3{\lambda }^2-{\beta }^2)\beta ch\left(\beta x\right)\mathit{sin}(\lambda x)+(3{\beta }^2-{\lambda }^2\left)\lambda sh\right(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)}}\hfill \\ b_{42}={\displaystyle \frac{(3{\lambda }^2-{\beta }^2)\beta ch\left(\beta x\right)\mathit{cos}(\lambda x)-(3{\beta }^2-{\lambda }^2\left)\lambda sh\right(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)}}\hfill \\ b_{13}={\displaystyle \frac{ch\left(\beta x\right)\mathit{sin}(\lambda x)}{2\beta \lambda EI}},b_{23}=-{\displaystyle \frac{ch\left(\beta x\right)\mathit{cos}(\lambda x)}{2\beta \lambda EI}},b_{33}=-{\displaystyle \frac{sh\left(\beta x\right)\mathit{sin}(\lambda x)}{2\beta \lambda EI}},b_{43}={\displaystyle \frac{sh\left(\beta x\right)\mathit{cos}(\lambda x)}{2\beta \lambda EI}}\hfill \\ b_{14}={\displaystyle \frac{\lambda ch\left(\beta x\right)\mathit{cos}(\lambda x)-\beta sh(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)EI}}, b_{24}={\displaystyle \frac{\lambda ch\left(\beta x\right)\mathit{sin}(\lambda x)+\beta sh(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)EI}}\hfill \\ b_{34}=-{\displaystyle \frac{\lambda sh\left(\beta x\right)\mathit{cos}(\lambda x)-\beta ch(\beta x\left)\mathit{sin}(\lambda x\right)}{2\beta \lambda \left({\lambda }^2+{\beta }^2\right)EI}}, b_{44}=-{\displaystyle \frac{\lambda sh\left(\beta x\right)\mathit{sin}(\lambda x)+\beta ch(\beta x\left)\mathit{cos}(\lambda x\right)}{2\beta \lambda ({\lambda }^2+{\beta }^2)EI}}\hfill \end{array}$$

(2)

When $N=2\sqrt{EIk}$

$$K_{xl}=\left[\begin{array}{cccc}x\mathit{sin}(\varphi x)& x\mathit{cos}(\varphi x)& \mathit{sin}(\varphi x)& \mathit{cos}(\varphi x)\\ \mathit{sin}(\varphi x)+\varphi x\mathit{cos}(\varphi x)& \mathit{cos}(\varphi x)-\varphi x\mathit{sin}(\varphi x)& \varphi \mathit{cos}(\varphi x)& -\varphi \mathit{sin}(\varphi x)\\ EI\varphi \left\{\varphi x\mathit{sin}(\varphi x\right)-2\mathit{cos}(\varphi x\left)\right\}& EI\varphi \left\{\varphi x\mathit{cos}(\varphi x\right)+2\mathit{sin}(\varphi x\left)\right\}& EI{\varphi }^2\mathit{sin}(\varphi x)& EI{\varphi }^2\mathit{cos}(\varphi x)\\ EI{\varphi }^2\left\{3\mathit{sin}(\varphi x\right)+\varphi x\mathit{cos}(\varphi x\left)\right\}& EI{\varphi }^2\left\{3\mathit{cos}(\varphi x\right)-\varphi x\mathit{sin}(\varphi x\left)\right\}& EI{\varphi }^3\mathit{cos}(\varphi x)& -EI{\varphi }^3\mathit{sin}(\varphi x)\end{array}\right]$$

and its inverse matrix is as follows:

$$K_{xl}^{-1}=\left[\begin{array}{cccc}{\displaystyle \frac{\varphi \mathit{cos}(\varphi x)}{2}}& -{\displaystyle \frac{\mathit{sin}(\varphi x)}{2}}& -{\displaystyle \frac{\mathit{cos}(\varphi x)}{2\varphi EI}}& {\displaystyle \frac{\mathit{sin}(\varphi x)}{2{\varphi }^{2}EI}}\\ -{\displaystyle \frac{\varphi \mathit{sin}(\varphi x)}{2}}& -{\displaystyle \frac{\mathit{cos}(\varphi x)}{2}}& {\displaystyle \frac{\mathit{sin}(\varphi x)}{2\varphi EI}}& {\displaystyle \frac{\mathit{cos}(\varphi x)}{2{\varphi }^{2}EI}}\\ {\displaystyle \frac{2\mathit{sin}(\varphi x)-\varphi x\mathit{cos}(\varphi x)}{2}}& {\displaystyle \frac{3\mathit{cos}(\varphi x)+\varphi x\mathit{sin}(\varphi x)}{2\varphi }}& {\displaystyle \frac{x\mathit{cos}(\varphi x)}{2\varphi EI}}& -{\displaystyle \frac{\mathit{cos}(\varphi x)+\varphi x\mathit{sin}(\varphi x)}{2{\varphi }^{3}EI}}\\ {\displaystyle \frac{2\mathit{cos}(\varphi x)+\varphi x\mathit{sin}(\varphi x)}{2}}& -{\displaystyle \frac{3\mathit{sin}(\varphi x)-\varphi x\mathit{cos}(\varphi x)}{2\varphi }}& -{\displaystyle \frac{x\mathit{sin}(\varphi x)}{2\varphi EI}}& {\displaystyle \frac{\mathit{sin}(\varphi x)-\varphi x\mathit{cos}(\varphi x)}{2{\varphi }^{3}EI}}\end{array}\right]$$

(3)

When $N>2\sqrt{EIk}$

$$K_{xl}=\left[\begin{array}{cccc}\mathit{sin}(\varphi x)& \mathit{cos}(\varphi x)& \mathit{sin}(\psi x)& \mathit{cos}(\psi x)\\ \varphi \mathit{cos}(\varphi x)& -\varphi \mathit{sin}(\varphi x)& \psi \mathit{cos}(\psi x)& -\psi \mathit{sin}(\psi x)\\ EI{\varphi }^2\mathit{sin}(\varphi x)& EI{\varphi }^2\mathit{cos}(\varphi x)& EI{\psi }^2\mathit{sin}(\psi x)& EI{\psi }^2\mathit{cos}(\psi x)\\ EI{\varphi }^3\mathit{cos}(\varphi x)& -EI{\varphi }^3\mathit{sin}(\varphi x)& EI{\psi }^3\mathit{cos}(\psi x)& -EI{\psi }^3\mathit{sin}(\psi x)\end{array}\right]$$

and its inverse matrix is as follows:

$$K_{xl}^{-1}=\left[\begin{array}{cccc}-{\psi }^2\mathit{sin}(\varphi x)& -{\displaystyle \frac{{\psi }^2\mathit{cos}(\varphi x)}{\varphi }}& {\displaystyle \frac{\mathit{sin}(\varphi x)}{EI}}& {\displaystyle \frac{\mathit{cos}(\varphi x)}{\varphi EI}}\\ -{\psi }^2\mathit{cos}(\varphi x)& {\displaystyle \frac{{\psi }^2\mathit{sin}(\varphi x)}{\varphi }}& {\displaystyle \frac{\mathit{cos}(\varphi x)}{EI}}& -{\displaystyle \frac{\mathit{sin}(\varphi x)}{\varphi EI}}\\ {\varphi }^2\mathit{sin}(\psi x)& {\displaystyle \frac{{\varphi }^2\mathit{cos}(\psi x)}{\psi }}& -{\displaystyle \frac{\mathit{sin}(\psi x)}{EI}}& -{\displaystyle \frac{\mathit{cos}(\psi x)}{\psi EI}}\\ {\varphi }^2\mathit{cos}(\psi x)& -{\displaystyle \frac{{\varphi }^2\mathit{sin}(\psi x)}{\psi }}& -{\displaystyle \frac{\mathit{cos}(\psi x)}{EI}}& {\displaystyle \frac{\mathit{sin}(\psi x)}{\psi EI}}\end{array}\right]/({\varphi }^2-{\psi }^2)$$

4. Solutions to the Influence Line Equations of the Main Beam

4.1. Composition of the Vertical Deflection of the Main Beam

By comparing the right terms of the equilibrium differential equations of the hybrid section, the pure suspension section, and the cable-stayed section in Equations (6), (9) and (15), respectively, it is found that the hybrid and pure suspension sections are equivalent to a concentrated load p and a uniformly distributed load $H_p{y}^{″}$ acting on the beam supported on the elastic foundation, while the cable-stayed section is equivalent to a concentrated load p acting on the beam supported on the elastic foundation.

When a vertical load p acts at the position of $\xi$ of the main beam, as shown in Figure 5a, the hybrid bridge system can be simplified as a concentrated load acting on the elastic foundation and equivalent uniformly distributed load $H_p{y}^{″}$ acting within the pure suspension section and the hybrid section, as shown in Figure 5b.

For the beam on the elastic foundation, assuming that the vertical beam deflection is ${\nu }_1(x,\xi )$ and ${\nu }_2(x,\xi )$, under the concentrated load P and equivalent uniform load $H_p{y}^{″}$, respectively, then the vertical deflection of the beam of the hybrid system under the concentrated load is as follows:

$$\nu \left(x,\xi \right)={\nu }_1\left(x,\xi \right)+{\nu }_2\left(x,\xi \right)$$

(17)

Based on the general solution of the above differential equations for different sections, the vertical deflection of the beam on the elastic foundation under the concentrated load P can be computed based on the following:

$${\nu }_1\left(x,\xi \right)=A\left(\xi \right)f_1\left(x\right)+B\left(\xi \right)f_2\left(x\right)+C\left(\xi \right)f_3\left(x\right)+D\left(\xi \right)f_4\left(x\right)$$

(18)

4.2. Solution to the Beam Equation on Elastic Foundation Under the Concentrated Load

Based on the derivation of the differential equations for various sections of the hybrid bridge system, the coefficients of the differential equations are all assumed to be constant, i.e., the cross-sectional properties, bending stiffness, elastic support stiffness, and axial force in each section of the main beam are constant. In fact, the stiffness and axial force of the elastic support of the main beam are constantly changing at different positions, which makes the equilibrium differential equation become a homogeneous differential equation with variable coefficients, and the general solution of the equation cannot be obtained. In order to solve the differential equations with variable coefficients, the main beam is divided into several sections. In each section, the bending stiffness and axial force of the main beam and the stiffness of the elastic support of the cable are taken as constants. Joints should also be formed at the point where the concentrated load is applied, and the beam section is divided into left and right sections to ensure that there is no additional load in the main beam section, and the additional load is all applied to the joints.

As shown in Figure 6, the main beam is discretized into n sections, the left-end ordinate of the section is $x_i$, the right-end ordinate is $x_{i+1}$, and n + 1 nodes are formed on the main beam, and the ordinate is in the sequence $x_1,x_2,\cdots ,x_i,x_{i+1},\cdots ,x_n,x_{n+1}$. The eigenvalue matrix of the section is $K_i\left(x\right)$, and the inverse matrix is $K_i^{-1}\left(x\right)$.

There are two main methods to solve the equation: the Transfer Matrix Method [32] and the Large Matrix Method. The Transfer Matrix Method is used to establish the relationship between the coefficient vector of any beam segment and the coefficient vector of a particular beam segment via matrix transfer, thus reducing the number of unknowns. Its advantage is that the unknown quantity is less, while its disadvantage is that the element value of the eigenvector of each beam segment is quite different, which will lead to the accumulation of errors in the transmission process, thus affecting the calculation accuracy. The Large Matrix Method is similar to the matrix displacement method in structural mechanics. Its characteristic is that there is an unknown coefficient vector in each beam segment, and the relationship between the unknown coefficients between two adjacent beam segments is determined by the joint boundary conditions at the common nodes of the adjacent beam segments. Its advantage is that there is no need to transfer between eigenvectors and the error caused by data accuracy is small. The disadvantage is that there are relatively many unknowns and the size of the matrix is large. In this study, the Transfer Matrix Method is employed, which is described as follows.

It is assumed that a concentrated load is applied at the location $x_k$ and a vertical constraint is set at the location $x_m$. Then, the characteristic mechanical vectors at the common endpoint of the main beam are equal except for the two positions $x_m$ and $x_k$. According to the principle that the mechanical characteristic quantities of the adjacent beam segments of the joint $x_i$ are equal, the following equation is obtained:

$$K_{i-1}\left(x_i\right)S_{i-1}=K_i\left(x_i\right)S_i$$

(19)

The coefficient vector of the adjacent beam segments of the node $x_i$ is as follows:

$$S_{i-1}={K_{i-1}^{-1}\left(x_i\right)K}_i\left(x_i\right)S_i$$

(20)

or

$$S_i={K_i^{-1}\left(x_i\right)K}_{i-1}\left(x_i\right)S_{i-1}$$

(21)

For the joint $x_k$ with the vertical concentrated load, the equilibrium equation at the left and right sides $x_k$ of the beams is as follows:

$$K_{k-1}\left(x_k\right)S_{k-1}=K_k\left(x_k\right)S_k+F_k$$

(22)

where $F_k=(0 0 0 P{)}^T$ is the additional load vector at this position. Then, the coefficient vector of the adjacent beam segment of the node $x_k$ is: as follows

$$S_{k-1}={K_{k-1}^{-1}\left(x_k\right)K}_k\left(x_k\right)S_k+K_{k-1}^{-1}\left(x_k\right)F_k$$

(23)

or

$$S_k={K_k^{-1}\left(x_k\right)K}_{k-1}\left(x_k\right)S_{k-1}-K_k^{-1}\left(x_k\right)F_k$$

(24)

For the joint $x_m$ with the vertical constraint, the equilibrium equation at the left and right sides $x_m$ of the beams is as follows:

$$K_{m-1}\left(x_m\right)S_{m-1}=K_m\left(x_m\right)S_m+F_m$$

(25)

where $F_m=(0 0 0 R_m{)}^T$ is the additional load vector at the location of $x_m$ and $R_m$ is the vertical force generated by the vertical constraint on the main beam, which is an unknown quantity. Then, the coefficient vector of the adjacent beam segment of the node $x_m$ is as follows:

$$S_{m-1}={K_{m-1}^{-1}\left(x_m\right)K}_m\left(x_m\right)S_m+K_{m-1}^{-1}\left(x_m\right)F_m$$

(26)

or

$$S_m={K_m^{-1}\left(x_m\right)K}_{m-1}\left(x_m\right)S_{m-1}-K_m^{-1}\left(x_m\right)F_m$$

(27)

For the Equations (20), (23) and (26), the relationship between the coefficient vector of the beam segment t and the coefficient vector of the beam segment n is as follows:

$$S_t=\left\{\begin{array}{ccc}{C}_{t,n}{S}_n\hfill & & k\le t<n\\ C_{t,n}{S}_n+C_{t,k-1}{K}_{k-1}^{-1}\left(x_k\right)F_k\hfill & & m\le t<k\\ C_{t,n}{S}_n+C_{t,k-1}{K}_{k-1}^{-1}\left(x_k\right)F_k+C_{t,m-1}{K}_{m-1}^{-1}\left(x_m\right)F_m\hfill & & 1\le t<m\end{array}\right.$$

(28)

where $C_{t,n}=\prod _{i=t+1}^n\left\{{K_{i-1}^{-1}\left(x_i\right)K}_i\left(x_i\right)\right\}$, $C_{t,k-1}=\prod _{i=t+1}^{k-1}\left\{{K_{i-1}^{-1}\left(x_i\right)K}_i\left(x_i\right)\right\}$, $C_{t,m-1}=\prod _{i=t+1}^{m-1}\left\{{K_{i-1}^{-1}\left(x_i\right)K}_i\left(x_i\right)\right\}$, and let $C_{t,t}=E$.

Or from Equations (21), (24) and (27), the relationship between the coefficient vector of the beam segment t and the coefficient vector of the beam segment 1 is as follows:

$$S_t=\left\{\begin{array}{ccc}{D}_{t,1}{S}_1\hfill & & 1\le t<m\\ D_{t,1}{S}_1-D_{t,m}{K}_m^{-1}\left(x_m\right)F_m\hfill & & m\le t<k\\ D_{t,1}{S}_1-D_{t,m}{K}_m^{-1}\left(x_m\right)F_m-D_{t,k}{K}_k^{-1}\left(x_k\right)F_k\hfill & & k\le t\le n\end{array}\right.$$

(29)

where $D_{t,1}=\prod _{i=t}^2\left\{{K_i^{-1}\left(x_i\right)K}_{i-1}\left(x_i\right)\right\}$,$D_{t,m}=\prod _{i=t}^{m+1}\left\{{K_i^{-1}\left(x_i\right)K}_{i-1}\left(x_i\right)\right\}$, $D_{t,k}=\prod _{i=t}^{k+1}\left\{{K_i^{-1}\left(x_i\right)K}_{i-1}\left(x_i\right)\right\}$, and let $D_{t,t}=E$.

According to the constraint equations of the left and right ends of the main beam, four linear equations with five unknowns can be formed. Due to the vertical constraint at the node $x_m$, an unknown quantity $R_m$ is added, and an independent constraint equation is added. The vertical displacement of the first element in the vector $K_m\left(x_m\right)S_m$ is 0. Five unknowns correspond to five equations. By solving the linear equations, the coefficients of each segment and the vertical force at the middle support can be obtained. Considering the typical situation that the left and right ends of the main beam are only set as the vertical constraint as an example, then the characteristic mechanical equations of the left and right ends of the main beam are as Equations (30) and (31), respectively.

$${\left[0 {v\left(x_1\right)}^{\prime } 0 Q\left(x_1\right)\right]}^T=K_1\left(x_1\right)S_1$$

(30)

$${\left[0 {v\left(x_{n+1}\right)}^{\prime } 0 Q\left(x_{n+1}\right)\right]}^T=K_n\left(x_{n+1}\right)({D_{n,1}S}_1-D_{n,m}{K}_m^{-1}\left(x_m\right)F_m-D_{n,k}{K}_k^{-1}\left(x_k\right)F_k)$$

(31)

The characteristic mechanical equation of the right end of $x_m$ for the main beam is as follows:

$${\left[0 {v\left(x_m\right)}^{\prime } M\left(x_m\right) Q\left(x_m\right)\right]}^T=K_m\left(x_m\right)({D_{m,1}S}_1-D_{m,m}{K}_m^{-1}\left(x_m\right)F_m)$$

(32)

Take the first and third rows of Equations (30) and (31) and the first row of Equation (32), it can form five linear equations with five unknowns, which consist of the coefficients of $A_1$, $B_1$, $C_1$, and $D_1$ of the vector $S_1$, and $R_m$ of the vector $F_m$. Then, the unique solution can be solved.

When the bending moment is acting on the main beam node k along with the vertical force, then the force vector is $F_k={\left[0 0 M P\right]}^T$.

From the above derivation, when the beam is supported at various locations, only the number of terms $x_m$ needs to be increased, and the number of unknown quantities and the equations of the final linear equations increases correspondingly according to the number of supports. When there are multiple known loads in the middle of the main beam, the corresponding terms $F_k$ need to be added, but the unknown quantity and the number of equations of the final linear equation system remain unchanged.

4.3. Solution to the Beam Equation Based on Elastic Foundation Under the Equivalent Uniformly Distributed Load

In order to solve the deflection ${\nu }_2\left(x,\xi \right)$ of the beam on the elastic foundation under the equivalent uniformly-distributed load $H_p{y}^{″}$, the horizontal force component of the main cable $H_p$ under the concentrated load p should be determined at first, and ${\nu }_2\left(x,\xi \right)$ could be obtained based on the integral form of the concentrated load.

$$v_2(x,\xi )=H_p{y}^{″}\int v_1(x,t)dt=H_p(y_0^{″}{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1^{″}{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3^{″}{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1(x,t)dt$$

(33)

For the purpose of simplification, the equation ${\int }_{x_1}^{x_2}f\left(x\right)dx+{\int }_{x_3}^{x_4}f\left(x\right)dx=({\int }_{x_1}^{x_2}+{\int }_{x_3}^{x_4}\left)f\right(x)dx$ is used.

It is seen that the Equation (33) consists of two unknowns, ${\nu }_2\left(x,\xi \right)$ and $H_p$. Therefore, the compatibility equation of the main cable should be used jointly.

$$\sum _{i=1}^n({\displaystyle \frac{H_p{L}_{\theta }}{E_c{A}_c}}+a_t{\Delta }_T{L}_t-\int y{\left(x\right)}^{\prime }w{(x,\xi )}^{\prime }dx)=0$$

(34)

where $L_{\theta }=\int {\displaystyle \frac{1}{{cos}^3\theta }}dx$ and $L_t=\int {\displaystyle \frac{1}{{cos}^2\theta }}dx$. $E_c$, $A_c$, $a_t$, ${\Delta }_t$, and $w\left(x,\xi \right)$ are the elastic modulus, cross-sectional area, coefficient of thermal expansion, temperature variation, and vertical deflection of the main cable, respectively.

4.3.1. Vertical Deflection Curve of the Main Cable in the Main Span

The main cable in middle span is divided into three sections, corresponding to the hybrid section, pure suspension section, and cable-stayed section, as shown in Figure 5a. The no-hanger main cable only bears its own uniform dead load $q_2$; therefore, it has $y_2^{″}=-{\displaystyle \frac{q_2}{H_q}}$, $y_2^{″}+w^{″}=-{\displaystyle \frac{q_2}{H_q+H_p}}$, where $y_2$ is the left no-hanger main cable’s initial shape. By subtracting these two equations, the vertical deflection of the right end of no-hanger main cable caused by the load is $w^{″}=-{\displaystyle \frac{H_p{y}_2^{″}}{H_q+H_p}}\approx {\displaystyle \frac{H_p{y}_2^{″}}{H_q}}$. For the right no-hanger main cable, $y_2^{″}=y_4^{″}$, where $y_4$ is the right no-hanger main cable’s shape function under the dead load. Similarly, the vertical deflection of the left end of the no-hanger main cable caused by the load is $w^{″}=-{\displaystyle \frac{H_p{y}_4^{″}}{H_q+H_p}}\approx {\displaystyle \frac{H_p{y}_4^{″}}{H_q}}={\displaystyle \frac{H_p{y}_2^{″}}{H_q}}$.

Therefore, the deflection curve of the main cable $w\left(x,\xi \right)$ is as follows:

$$w(x,\xi )=\left\{\begin{array}{cc}{\displaystyle \frac{H_p{y}_2^{″}}{2H_q}}({\displaystyle \frac{L}{2}}+x)(L_{xl}-{\displaystyle \frac{L}{2}}-x)+{\displaystyle \frac{v(L_{xl}-\frac{L}{2},\xi )}{L_{xl}}}({\displaystyle \frac{L}{2}}+x)& -{\displaystyle \frac{L}{2}}\le x<L_{xl}-{\displaystyle \frac{L}{2}}\\ v(x,\xi )& L_{xl}-{\displaystyle \frac{L}{2}}\le x\le {\displaystyle \frac{L}{2}}-L_{xl}\\ {\displaystyle \frac{H_p{y}_2^{″}}{2H_q}}(x-{\displaystyle \frac{L}{2}}+L_{xl})({\displaystyle \frac{L}{2}}-x)+{\displaystyle \frac{v(\frac{L}{2}-L_{xl},\xi )}{L_{xl}}}({\displaystyle \frac{L}{2}}-x)& {\displaystyle \frac{L}{2}}-L_{xl}<x\le {\displaystyle \frac{L}{2}}\end{array}\right.$$

(35)

4.3.2. Solution of the Compatibility Equation

For the hybrid bridge, the longitudinal constraints on the top of the main tower, including the main tower, stay cables, and main cables of the side span, can be simulated with an equivalent longitudinal stiffness $K_{tx}$ at the two main cable saddles, as shown in Figure 7.

When the equivalent longitudinal stiffness $K_{tx}$ is used, the Equation (34) can be rewritten as follows:

$${\displaystyle \frac{H_p}{K_{tx,1}}}+{\displaystyle \frac{H_p}{K_{tx,2}}}+{\displaystyle \frac{H_p{L}_{\theta }}{E_c{A}_c}}+a_t{\Delta }_T{L}_t-{\int }_{-{\displaystyle \frac{L}{2}}}^{{\displaystyle \frac{L}{2}}}y{\left(x\right)}^{\prime }w{(x,\xi )}^{\prime }dx=0$$

(36)

Integrating the last term of the left side of Equation (36) and consider the initial conditions $w\left(-{\displaystyle \frac{L}{2}}, \xi \right)=0$ and $w\left({\displaystyle \frac{L}{2}}, \xi \right)=0$, it gives the following:

$${\int }_{-{\displaystyle \frac{L}{2}}}^{{\displaystyle \frac{L}{2}}}y{(x)}^{\prime }w{(x,\xi )}^{\prime }dx=[y(x)\prime w(x,\xi )]{|}_{x=-{\displaystyle \frac{L}{2}}}^{x={\displaystyle \frac{L}{2}}}-\sum _{i=0}^4{y}_i^{″}\int w(x,\xi )dx=-\sum _{i=0}^4{y}_i^{″}\int w(x,\xi )dx$$

(37)

where

$$\begin{array}{c}{y}_0^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}w(x,\xi )dx}=y_0^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{v}_1(x,\xi )dx}+H_p(y_0^{″2}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_0^{″}{y}_3^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}{y}_0^{″}{y}_1^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}})v_1(x,t)dt}}dx\hfill \\ y_1^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}w(x,\xi )dx}=y_1^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{v}_1(x,\xi )dx}+H_p(y_1^{″}{y}_0^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_1^{″}{y}_3^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}{y}_1^{″2}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}})v_1(x,t)dt}}dx\hfill \\ y_3^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}w(x,\xi )dx}=y_3^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{v}_1(x,\xi )dx}+H_p(y_3^{″}{y}_0^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_3^{″2}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}{y}_3^{″}{y}_1^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}})v_1(x,t)dt}}dx\hfill \\ y_2^{″}{\int }_{{\displaystyle \frac{L}{2}}-L_{xl}}^{{\displaystyle \frac{L}{2}}}w(x,\xi )dx={\displaystyle \frac{H_p{y_2^{″}}^2{L}_{xl}^3}{12H_q}}+{\displaystyle \frac{y_2^{″}{L}_{xl}}{2}}v({\displaystyle \frac{L}{2}}-L_{xl},\xi )\hfill \\ v({\displaystyle \frac{L}{2}}-L_{xl},\xi )=v_1({\displaystyle \frac{L}{2}}-L_{xl},\xi )+H_p(y_0^{″}{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1^{″}{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3^{″}{\int }_{-{\displaystyle \frac{L_{xd}}{2}}-L_{jc}}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1({\displaystyle \frac{L}{2}}-L_{xl},t)dt\hfill \\ y_4^{″}{\int }_{-{\displaystyle \frac{L}{2}}}^{L_{xl}-{\displaystyle \frac{L}{2}}}w(x,\xi )dx={\displaystyle \frac{H_p{y_2^{″}}^2{L}_{xl}^3}{12H_q}}+{\displaystyle \frac{y_2^{″}{L}_{xl}}{2}}v(L_{xl}-{\displaystyle \frac{L}{2}},\xi )\hfill \\ v(L_{xl}-{\displaystyle \frac{L}{2}},\xi )=v_1(L_{xl}-{\displaystyle \frac{L}{2}},\xi )+H_p(y_0^{″}{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1^{″}{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3^{″}{\int }_{-{\displaystyle \frac{L_{xd}}{2}}-L_{jc}}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1(L_{xl}-{\displaystyle \frac{L}{2}},t)dt\hfill \end{array}$$

Rewrite the Equation (37) and let

$$\begin{gathered} {\Lambda }_1=(y_0^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}{y}_1^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}+}{y}_3^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}})v_1(x,\xi )dx}+{\displaystyle \frac{y_2^{″}{L}_{xl}}{2}}{v}_1({\displaystyle \frac{L}{2}}-L_{xl})+{\displaystyle \frac{y_2^{″}{L}_{xl}}{2}}{v}_1(L_{xl}-{\displaystyle \frac{L}{2}}) \\ \begin{array}{ll}{\Lambda }_2& =(y_0^{″2}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_0^{″}{y}_3^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}{y}_0^{″}{y}_1^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}+}}{y}_1^{″}{y}_0^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_1^{″}{y}_3^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}\\ & y_1^{″2}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}+}}{y}_3^{″}{y}_0^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}}{y}_3^{″2}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}+}}{y}_3^{″}{y}_1^{″}{\displaystyle {\int }_{-(\frac{L_{xd}}{2}+L_{jc})}^{-\frac{L_{xd}}{2}}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}})v_1(x,t)dt}}dx+\\ & {\displaystyle \frac{y_2^{″2}{L}_{xl}^3}{6H_q}}+{\displaystyle \frac{y_2^{″}{L}_{xl}}{2}}(y_0^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}}+}{y}_1^{″}{\displaystyle {\int }_{\frac{L_{xd}}{2}}^{\frac{L_{xd}}{2}+L_{jc}}+}{y}_3^{″}{\displaystyle {\int }_{-\frac{L_{xd}}{2}-L_{jc}}^{-\frac{L_{xd}}{2}})[v_1(\frac{L}{2}-L_{xl},t)+v_1(L_{xl}-\frac{L}{2},t)]dt},\end{array} \end{gathered}$$

This gives the following:

$$\sum _{i=0}^4{y}_i^{″}\int w(x,\xi )dx={\Lambda }_1+{\Lambda }_2{H}_p$$

(38)

Substitute Equation (38) into Equation (36), then the horizontal force component of the main cable is solved.

$$H_p=-{\displaystyle \frac{a_t{\Delta }_T{L}_t+{\Lambda }_1}{\frac{1}{K_{tx,1}}+\frac{1}{K_{tx,2}}+\frac{L_{\theta }}{E_c{A}_c}+{\Lambda }_2}}$$

(39)

Substitute Equation (39) into Equation (33), the vertical deflection ${\nu }_2\left(x,\xi \right)$ of the beam on the elastic foundation under the equivalent uniformly distributed load can be solved as follows:

$$v_2(x,\xi )={\displaystyle \frac{a_t{\Delta }_T{L}_t+{\Lambda }_1}{\frac{1}{K_{tx,1}}+\frac{1}{K_{tx,2}}+\frac{L_{\theta }}{E_c{A}_c}+{\Lambda }_2}}(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1(x,t)dt$$

(40)

4.4. Influence Lines of the Main Beam

The basic equation of the deflection ${\nu }_1\left(x,\xi \right)$ of the beam on the elastic foundation under the concentrated load is $v_1(x,\xi )=A\left(\xi \right)f_1\left(x\right)+B\left(\xi \right)f_2\left(x\right)+C\left(\xi \right)f_3\left(x\right)+D\left(\xi \right)f_4\left(x\right)$.

Under the action of the concentrated load, by substituting Equation (33), the deflection curve $\nu \left(x,\xi \right)={\nu }_1\left(x,\xi \right)+{\nu }_2\left(x,\xi \right)$ of the main beam of the hybrid system is as follows:

$$v(x,\xi )=v_1(x,\xi )+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1(x,t)dt$$

(41)

Among them, the integral of deflection ${\nu }_1\left(x,t\right)$ relative to the load position $t$ can be written in the following form:

$$\int v_1(x,t)dt=\int A\left(t\right)dt{f}_1\left(x\right)+\int B\left(t\right)dt{f}_2\left(x\right)+\int C\left(t\right)dt{f}_3\left(x\right)+\int D\left(t\right)dt{f}_4\left(x\right)$$

(42)

Double integrating the ${\nu }_1\left(x,t\right)$ relative to the load position and longitudinal coordinate, it gives the following:

$$\iint v_1(x,t)dtdx=\int A\left(t\right)dt\int f_1\left(x\right)dx+\int B\left(t\right)dt\int f_2\left(x\right)dx+\int C\left(t\right)dt\int f_3\left(x\right)dx+\int D\left(t\right)dt\int f_4\left(x\right)dx$$

(43)

Based on the Equation (41), the derivative of the deflection of the main beam, the moment, and shear force of the main beam for the hybrid system under the concentrated load can be computed as follows:

$$v{(x,\xi )}^{\prime }=v_1{(x,\xi )}^{\prime }+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})v_1{(x,t)}^{\prime }dt$$

(44)

$$M(x,\xi )=M_1(x,\xi )+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})M_1(x,t)dt$$

(45)

$$Q(x,\xi )=Q_1(x,\xi )+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})Q_1(x,t)dt$$

(46)

where ${\nu }_1\left(x,\xi \right)$, ${\nu }_1{\left(x,\xi \right)}^{\prime }$, $M_1\left(x,\xi \right)$, and $Q_1\left(x,\xi \right)$ are the deflection, derivative of deflection, moment, and shear force of the beam supported on the elastic foundation. They form the mechanical characteristics vector $T_1=({\nu }_1\left(x,\xi \right)$, ${\nu }_1{\left(x,\xi \right)}^{\prime }$, $M_1\left(x,\xi \right)$, $Q_1\left(x,\xi \right){)}^T$. Based on Section 4.2, the coefficient vector $S=(A B C D{)}^T$ can be solved. Then, the mechanical characteristic vector T can be solved as $T=KS$.

Equations (41)–(46) can be rewritten as the matrix form.

$$T=K[S+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})Sdt]$$

(47)

Or

$$\left\{\begin{array}{c}v(x,\xi )\\ v{(x,\xi )}^{\prime }\\ M(x,\xi )\\ Q(x,\xi )\end{array}\right\}=K(x)[\left\{\begin{array}{c}A(\xi )\\ B(\xi )\\ C(\xi )\\ D(\xi )\end{array}\right\}+H_p(y_0″{\int }_{-{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}}+y_1″{\int }_{{\displaystyle \frac{L_{xd}}{2}}}^{{\displaystyle \frac{L_{xd}}{2}}+L_{jc}}+y_3″{\int }_{-({\displaystyle \frac{L_{xd}}{2}}+L_{jc})}^{-{\displaystyle \frac{L_{xd}}{2}}})\left\{\begin{array}{c}A(t)\\ B(t)\\ C(t)\\ D(t)\end{array}\right\}dt]$$

(48)

5. Case Study for Model Verification

5.1. Model Verification

In this section, a case study is presented to verify the proposed model by comparing it with the 2D numerical simulation results from the commercial software MIDAS Civil 2021. The case bridge is simplified from the Xihoumen Bridge, which is a combined highway and railway bridge with a main span length of 1488 m, and the shape of the main cable is catenary. The schematic of the hybrid bridge system in Figure 5 is employed in this section. For the FEM model, the initial stage is determined based on the completion stage of the bridge under the permanent dead load to calculate the stress-free stress and coordinates of the main cable and hangers. As shown in Figure 5a, the parameters of the hybrid bridge are as follows. The lengths of the pure suspension section $L_{xd}$, the hybrid section at each side $L_{jc}$, the cable-stayed section at each side $L_{xl}$, the auxiliary span $L_f$, and the side span $L_b$ are 364 m, 168 m, 350 m, 406 m, and 112 m, respectively. The rise-to-span ratio of the main cable is 1/6.5. Two main cables are used. The equivalent weight density of the main cable is ${\gamma }_c=82 \mathrm{k}\mathrm{N}/{\mathrm{m}}^3$. The total area of two main cables is $A_c=0.866 {\mathrm{m}}^2$, and hence, the main cable‘s own dead load is $q_c=71 \mathrm{k}\mathrm{N}/\mathrm{m}$, and the elastic modulus of the main cables is $E_c={2.0\times 10}^5 \mathrm{M}\mathrm{P}\mathrm{a}$. For the pure suspension section, the load density bear by hangers is $q_{h0}=577 \mathrm{k}\mathrm{N}/\mathrm{m}$, and the corresponding total load density of the main cables is $q_0=q_c+q_{h0}=641 \mathrm{k}\mathrm{N}/\mathrm{m}$. For the hybrid section, the load density taken by hangers is $q_{h1}=415 \mathrm{k}\mathrm{N}/\mathrm{m}$, and the corresponding total load density of the main cables is $q_1=q_c+q_{h1}=486 \mathrm{k}\mathrm{N}/\mathrm{m}$. Then, based on the parabolic approximation for the main cables, the horizontal load component of the main cables due to the dead load is $H_q=$ 518,175 kN.

The bending stiffness of the main beam is $EI=2.1\times {10}^{9 }\mathrm{k}\mathrm{N}/\mathrm{m}$. For the cable-stayed section and the auxiliary span and the side span, the elastic stiffness of the main beam provided by the cables is $k=400 \mathrm{k}\mathrm{N}/{\mathrm{m}}^2$. In the hybrid section, the elastic stiffness of the main beam provided by the cables is $k=200 \mathrm{k}\mathrm{N}/{\mathrm{m}}^2$. The axial load sustained by the main beam N is $1.5\times {10}^5 \mathrm{k}\mathrm{N}$, $4.0\times {10}^5 \mathrm{k}\mathrm{N}$, $4.0\times {10}^5 \mathrm{k}\mathrm{N}$, and $1.5\times {10}^{5 }\mathrm{k}\mathrm{N}$ in the side span, auxiliary span, cable-stayed section, and hybrid section, respectively. In the pure suspension section, the axial force of the main beam is neglected. The longitudinal constraint stiffness of the theoretical top of the main cables is $K_{tx}=2.5\times {10}^5 \mathrm{k}\mathrm{N}/\mathrm{m}$.

Considering a unit load of 1 kN applied in the middle-span, the deflection and moment of the elastic main beam are computed using both the numerical model and the proposed model. Comparisons of the deflection curve and moment are shown in Figure 8a and 8b, respectively. As seen in Figure 8a, the deflection curve computed using the proposed model matches well with that by the numerical model. The maximum deflection of the main beam at the middle-span is $1.53\times {10}^{-4} \mathrm{m}$ and $1.45\times {10}^{-4} \mathrm{m}$ based on the numerical model and the proposed model, respectively, and the difference is −5%, which is small. Similarly, as seen in Figure 8b, the moment computed using the proposed model also matches well with that using the numerical model. The maximum moment of the main beam at the middle-span is $30.86 \mathrm{k}\mathrm{N}·\mathrm{m}$ and $30.73 \mathrm{k}\mathrm{N}·\mathrm{m}$ based on the numerical model and the proposed model, respectively, and the difference is −0.4%, which is also small. Therefore, when a concentrated unit load is applied at the middle span, the proposed method can give an accurate estimation of the deflection and moment of the main beam.

5.2. Characteristics of Influence Lines

When a unit load 1 kN is applied at an arbitrary location $\xi \in [-1218, 1218] \mathrm{m}$, the corresponding vertical deflection, rotation, moment, and shear of the main beam at a point $x\in [-1218, 1218] \mathrm{m}$ is computed using the proposed model. Then, the influence lines of vertical deflection, rotation, moment, and shear of the main beam at all $x\in [-1218, 1218] \mathrm{m}$ are computed and the results are shown in Figure 9a, 9b, 9c, and 9d, respectively.

As seen in Figure 9a, the larger values of the vertical deflection influence line are symmetrical about the vertical axis at $\xi =x$. When the unit load is applied at $\xi =0 \mathrm{m}$, the maximum value of the vertical deflection influence line is $8.01\times {10}^{-5} \mathrm{m},$ which occurs for the influence line of the middle-span $x=0$. Meanwhile, the values of the vertical deflection influence lines within the pure suspension section are obviously larger than the hybrid and cable-stayed sections.

Similarly, in Figure 9b, the larger values of the beam rotation influence the line distributed at the two sides of $\xi =x$, and the distribution range is larger than that of the vertical deflection. When the unit load is applied at $\xi =-84 \mathrm{m}$, the maximum value of the beam rotation influence line is $4.98\times {10}^{-7} rad$, which occurs for the influence line at $x=-161 \mathrm{m}$.

It is seen in Figure 9c that the beam moment influence lines are all approximately symmetric about $\xi =x$, and its distribution range is larger than that of the vertical deflection. When the unit load is applied at $\xi =35 \mathrm{m}$, the maximum value of the beam rotation influence line is $25.35 \mathrm{k}\mathrm{N}·\mathrm{m}$ which occurs at the concerned location $x=35 \mathrm{m}$.

As for the beam shear in Figure 9d, the influence lines are symmetric about $\xi =x$. When the load is applied around the vertical constraint, the value of the beam shear influence line is relatively larger. When the unit load is applied at the auxiliary pier $\xi =-1106 \mathrm{m}$, the maximum value of the beam shear influence line is $0.98 \mathrm{k}\mathrm{N},$ which occurs for the influence line at $x=-1106 \mathrm{m}$.

Based on Figure 9, the influence line at an arbitrary location x can be easily obtained by cutting a cross profile at x. For a comparison purpose, the influence lines at six locations are obtained and plotted in Figure 10a, 10b, 10c, and 10d for the vertical deflection, rotation, moment, and shear of the main beam, respectively, which are at the middle-span $\mathrm{x}=0$, 1/4 of the pure suspension section $\mathrm{x}=91 \mathrm{m}$, the intersection between the suspension and hybrid sections $\mathrm{x}=182 \mathrm{m}$, the middle of the hybrid section $\mathrm{x}=266 \mathrm{m}$, the intersection between the hybrid and cable-stayed sections $\mathrm{x}=350 \mathrm{m}$, and the middle of the cable-stayed section $\mathrm{x}=525 \mathrm{m}$.

It is seen in Figure 10a that the vertical deflection influence line of the beam at any location is approximately symmetrical about the concerned location x, with the peak values occurring at the concerned location. In addition, affected by the deflection of the main cables, the vertical deflection influence line of the beam at any location has negative values at the two sides of the concerned point.

As seen in Figure 10b, the beam rotation influence line at any location is approximately anti-symmetrical with the concerned location. In addition, the peak values of the rotation influence lines at the pure suspension section are larger than those at the hybrid and cable-stayed sections.

For the beam moment in Figure 10c, the influence line at any concerned location, the line is approximately symmetrical about the concerned location, and the maximum value occurs at the concerned location. Moreover, the maximum values of the beam moment influence lines at the pure suspension section are larger than those at the hybrid and cable-stayed sections. However, the minimum values of the beam moment influence lines and its distribution range in the hybrid section are larger than those in the suspension and cable-stayed sections.

Similarly to the beam rotation in Figure 10b, the beam shear influence line at any location in Figure 10d is approximately anti-symmetrical about the concerned location. In addition, the maximum value of the shear influence line occurs at the concerned location. Moreover, the maximum values of the beam shear influence lines are close to the minimum values.

6. Conclusions

This paper presents a novel analytical model to analyze the mechanical performance of the cable-stayed suspension hybrid bridge, which is developed based on the modified suspension–elastic foundation beam theory with the consideration of the overlapped section between the suspension and stayed cables. The governing equations of the hybrid bridge system are derived based on the elastic beam foundation theory and the deflection theory, and the general solutions are also given. The influence lines of the main beam under the concentrated load are solved using the Transfer Matrix Method, while for the hybrid bridge system, the influence lines are solved by establishing the compatibility equation of the main cable with consideration of the longitudinal displacement of the main tower and no variation of unstressed length of the main cable. A case study is then presented to verify the proposed model using the finite element method. Meanwhile, the characteristics of the influence lines of the hybrid bridge system are investigated using the proposed model. The general conclusions are as follows.

(1)

The proposed model can accurately estimate the global deflection and internal forces of the hybrid bridge with the deviation within 5%. The derivation and solution of the proposed model can be easily complied in any programming platform with simple input parameters. Hence it can provide the theoretical basis for the verification of any complex finite element model with complex input data. In addition, the proposed model can overcome the difficulties in the finite element method during the bridge optimization, and hence, it can be used in the preliminary design stage.

(2)

For the complex hybrid system, the vertical deflection influence line of the beam at any location is approximately symmetrical about the concerned location, with the peak values occurring at the concerned location. In addition, a negative deflection zone for the main beam exists at a distance from the concentrated vertical load, which is mainly caused by the deflection of the main cables, leading to cambering of the beam.

(3)

The influence lines for the moment and vertical deflection of the beam are approximately symmetrical about the concerned locations, while those for beam rotation and shear are approximately anti-symmetric about the concerned locations. In addition, the maximum values of the influence lines for the vertical deflection and beam rotation at the pure suspension section are larger than those at the hybrid and cable-stayed sections.

It is clear that the proposed analytical model can accurately compute the structural deformation and internal forces of the long-span hybrid cable-stayed suspension bridge. However, due to some simplified assumptions made in the analytical model, such as the use of linear deflection theory and neglection of $H_p{\nu }^{″}$ in Equation (5), the results will, to some extent, deviate from the actual ones. In addition, due to the space limit, the methods to compute some critical parameters for the analytical model have not been given. These will be included in a companion paper to introduce the methods to compute the critical parameters, such as the longitudinal constraint stiffness of the main cable at the tower top, the vertical elastic support stiffness of the stayed cables to the main beam, and so on, and to perform a comprehensive parametric study on the structural performance of the bridge, including the span ratio of the cable-stayed part and sag ratio of the suspension part, and so on.