A Minimal GBT Model for Distortional-Twist Elastic Analysis of Box-Girder Bridges

Abstract

A simple and efficient method is proposed for the analysis of twist of rectangular box-girder bridges, which undergo distortion of the cross section. The model is developed in the framework of the Generalized Beam Theory and oriented towards semi-analytical solutions. Accordingly, only two modes are accounted for: (i) the torsional mode, in which the box-girder behaves as a Vlasov beam under nonuniform torsion, and, (ii) a distortional mode, in which the cross section behaves as a planar frame experiencing skew-symmetric displacements. By following a variational approach, two coupled, fourth-order differential equations in the modulating amplitudes are obtained. The order of magnitude of the different terms is analyzed, and further reduced models are proposed. A sample system, taken from the literature, is considered, for which generalized displacement and stress fields are evaluated. Both a Fourier solution for the coupled problem and a closed-form solution for the uncoupled problem are carried out, and the results are compared. Finally, the model is validated against finite element analyses.

1. Introduction

It is well-known that, when a mono-cellular box-girder undergoes torsion, induced by eccentric loads with respect to its longitudinal axis, its cross section suffers a distortion in its own plane, which modifies the original shape [1]. Such a phenomenon is often referred to, in the technical literature devoted to bridges, as differential flexure, since a quota of the external torsional moment not equilibrated by the internal torsional moment is instead bared by internal forces triggered by equal and opposite flexures of the two webs. Here, however, the word distortion is preferred. The phenomenon is dangerous, since flexure of the webs entails longitudinal normal stress, which adds to the tangential stresses of the Bredt theory of uniform torsion and to the normal stresses of the Vlasov theory of nonuniform torsion. Moreover, flanges and webs all suffer transverse bending, which would not be present if the cross section could maintain its shape. To limit such effects, diaphragms or bracings are occasionally introduced into the box-girder, but more often, these devices are omitted for construction simplicity. In these latter cases, a more refined analysis is required, in which the simple model of a beam must be abandoned and the box-girder must be modeled as a plate assembly.

The most popular approach to this problem is based on the so-called Beam on Elastic Foundation (BEF) analogy, attributed to Wright [2], after the pioneristic work by Vlasov. According to the BEF analogy, the effects of distortion are uncoupled from those of torsion and are governed by the classical differential equation of the beam on Winkler soil. The model was derived on a physical rather than a mathematical ground by substituting the continuous model of a beam with an assembly of elastic thin slices of the box-girder, rigidly connected among them. Among old papers, in Reference [3], an alternative physical model was proposed, leading to the same equation, and an example of reinforced concrete monocellular box-girder bridge was worked out; in Reference [4], parametric analyses were carried out and the results were given in different load conditions. The BEF analogy, in spite of its simplistic idea of uncoupling, is still used nowadays. As a few examples, in Reference [5], the analogy was exploited, aimed at supplying a useful tool for designers and accounting for the presence of diaphragms. In Reference [6], the approach is compared with finite element simulations, concluding that the simplified method provides a safe estimation of the response. In Reference [7], an extended version of the equation of BEF, including the second derivative, is used to analyze the influence of shear-deformable diaphragms.

In spite of using BEF, since 1990, it was clearly stated, in Reference [8], torsion and distortion cannot be separated since they are intrinsically coupled. The model presented there concerns a rectangular box-girder undergoing torsion–distorsion. It is based on simplified kinematics, in which three generalized configuration variables, depending on the longitudinal abscissa, are a priori introduced: one for twist, two (not independent) for bending of webs/flanges in their own planes, and one for distortion. The relevant strains are of membrane type, piece-wise linear, or constant along the directrix, while flexural and torsional plate contributions are ignored. By following an energy approach, a differential system of sixth-order in the configuration variable is derived and solved. In a successive paper [9], the author extended his theory to trapezoidal cross sections by introducing twelve configuration variables.

To analyze coupling effects, we need to resort to more refined models. Recently, higher-order theories of one-dimensional beams have been used, in which some added kinematic descriptors account for warping and distortion of the beam. An overview of different models in the literature is provided in Reference [10]. In this class of papers falls, e.g., References [11,12], in which warping and ovalization of a multi-layered beam is analyzed, and References [13,14,15], where nonlinear effects are accounted for.

The Finite Element Method (FEM) is the most popular tools today to solve problems. In the last decades, however, alternative numerical instruments have also been used, tailored to Thin-Walled Beam (TWB) structures, namely (i) the Finite Strip Method (FSM) [16,17,18] and (ii) the Generalized Beam Theory (GBT) [19]. The latter is a numerical version of the well-known Kantarovitch semi-variational method, in which the displacement field is taken as a linear combination of known functions defined on the cross-sectional domain (usually said modes) and unknown modulating amplitudes defined on the beam axis. The original two-dimensional plate problem is thus brought to a simpler one-dimensional problem in amplitude functions, for which Fourier series expansions, or 1D-FEM, can be applied. The idea had Vlasov himself as a forerunner, who used it to analyze closed TWB; some numerical applications of the Vlasov procedure have been recently illustrated in [20] by never mentioning GBT. The modern version of the method was implemented by many authors, following an original idea by Schardt, successively disseminated in English by References [19,21,22]. It had in Camotim and its school [23,24,25] a huge impulse in the last two decades, through a series of papers in which the choice of modes (cross-section analysis) and various algorithmic aspects were progressively refined. Luongo, Ranzi, and coworkers [26,27,28,29,30] launched the idea to combine GBT with a dynamic boundary value problem to automatize the cross-sectional analysis. De Miranda et al. [31,32] analyzed, in particular, the shear effects. More recently, in References [33,34], GBT was combined with FEM to solve problems in which localized strains and stresses occur. Moreover, the use of the Boundary Element Method (BEM) has also been proposed [35] to analyze beams with general (not necessarily thin-walled) beams.

All such philosophies, based on purely numerical approaches, are more oriented towards describing how a phenomenon manifests itself rather than towards explaining why it occurs, which are the essential parameters that govern them. Therefore, it is believed by the authors of this paper that ancient analytical methods, oriented towards hand-solving a problem, should be rethought in the light of the new knowledge and instruments available. In this view, a simplified model of the box-girder in torsion–distortion is proposed here. It is developed by reconsidering the approach of Reference [8] by framing it in the modern GBT. Accordingly, just two modes (and therefore two configuration variables) describing torsion and distortion, respectively, are shown to be accurate enough to capture the mechanical behavior based on an orthogonal decomposition of the planar cross-sectional displacement field. Strong differences from Reference [8], however, exist, namely (a) generalized configuration variables are not a priori assumed but, in the spirit of GBT, they naturally appear as the amplitude of the deformation modes; (b) distortion involves warping on a kinematic ground, entailing active normal stresses and passive tangential stresses (the reverse occurs in Reference [8]); (c) the plate behavior of the single elements is accounted for, allowing evaluation of the stresses within the thickness; (d) shear strains are not simply piece-wise constant on the directrix, allowing for satisfaction of the equilibrium at joints; and (e) an eighth-order (instead of sixth-order) differential system is derived, as a consequence of the accounted plate behavior and the different treatment of warping.

The main contributions of the paper are the following: (a) to provide a semi-analytical solution, oriented towards hand calculations, able to give reasonably accurate answers to a complex problem, which would otherwise require refined numerical models; (b) to perform an order of magnitude analysis of the different terms in the model, to formulate the simplest possible, and to gain insight into the relative importance of the distortion against the torsion in terms of geometrical parameters; (c) to explain that the effects of distortion can be described by just one mode, for which the in-plane component is taken orthogonal to that of torsion; (d) to clarify that coupling, in contrast, exists due to warping; and (e) to highlight that in an internally kinematically constrained model, the reactive stresses can play an important role. Overall, the paper is inspired by a methodological approach rather than by a computational philosophy. However, the ready-for-use formulas obtained here are believed to be of technical interest.

The paper is organized as follows. In Section 2, some background is supplied. In Section 3, the model is developed, and the governing equations are derived. In Section 4, some algorithmic aspects are addressed. In Section 5, a sample system is studied and the numerical results are discussed. In Section 6, some conclusions are drawn. A few Appendixes are devoted to illustrating details.

2. Background

Existing theories, necessary to formulate the model, are shortly summarized. They concern (i) the Vlasov theory [1], extended to closed TWB by Umanski (see, e.g., Reference [36]), and (ii) the Generalized Beam Theory [23,27], holding for plate assemblies. The box-girder refers to the longitudinal abscissa $z\in \left(0,\ell \right)$, running along the straight axis, and to the transverse abscissa s, running along the middle-line cross section $\Gamma$; the thickness $t\left(s\right)$, is assumed to be much smaller than the medium radius of the cross section.

2.1. Beam Theory

When a thin-walled box-girder is subject to nonuniform torsion, one can apply the extended Vlasov theory [36] for closed TWB. The theory assumes (i) that the cross section is undeformable in its own plane but warps out-of-plane and (ii) that the shear strains ${\gamma }_{zs}$ on the middle surface of the beam do not vanish, as for open TWB, but they satisfy the condition ${\gamma }_{zs}\left(z,s\right)t\left(s\right)=:\varphi \left(z\right)$, i.e., the flow of the strains across the thickness, at a given z, is constant along $\Gamma$, as suggested by the Bredt theory. The combination of the two geometric constraints yields the warping $w\left(z,s\right)=-{\theta }^{\prime }\left(z\right)\omega \left(s\right)$, with ${\theta }^{\prime }\left(z\right)$ being the genearally nonuniform torsional curvature and $\omega \left(s\right)$ being the warping function:

$$\omega \left(s\right)=2\left(\mathsf{\Omega }\left(s\right)-{\mathsf{\Omega }}_0\frac{{\int }_0^s\frac{ds}{t\left(s\right)}}{\oint \frac{ds}{t\left(s\right)}}\right)$$

(1)

where ${\mathsf{\Omega }}_0$ is the area enclosed by $\Gamma$ and $\mathsf{\Omega }\left(s\right)$ is the sector area function. Equation (1) generalizes the more familiar $\omega \left(s\right)=2\mathsf{\Omega }\left(s\right)$, which holds for open TWB.

When the equilibrium is expressed in terms of the twist angle $\theta \left(z\right)$, the following field equation and boundary conditions are found:

$$\begin{array}{cc}\hfill E{I}_{\omega }{\theta }^{″″}-GJ{\theta }^{″}& =c_z\left(z\right)\hfill \\ \hfill {\left[\left(GJ{\theta }^{\prime }-E{I}_{\omega }{\theta }^{‴}\right)\delta \theta \right]}_0^{\ell }& =0\hfill \\ \hfill {\left[E{I}_{\omega }{\theta }^{″}\delta {\theta }^{\prime }\right]}_0^{\ell }& =0\hfill \end{array}$$

(2)

where $E,G$ are elastic moduli; $I_{\omega }:={\oint }_{\Gamma }{\omega }^2\left(s\right)t\left(s\right)ds$ is the warping inertia moment, $J:=\frac{4{\mathsf{\Omega }}^2}{\oint \frac{ds}{t\left(s\right)}}$ is the De Saint Venant torsion inertia moment; $c_z\left(z\right)$ are distributed external torsional couples, $\delta$ indicates a virtual displacement, and a prime denotes z-differentiation. The second equation in Equation (2) shows that the external torsional moment ${\overline{M}}_t$ is equilibrated by two internal forces: the De Saint Venant torsional moment $M_t=GJ{\theta }^{\prime }$ and the complementary torsional moment $M_t^{\ast }=-E{I}_{\omega }{\theta }^{‴}$, i.e., ${\overline{M}}_t=M_t+M_t^{\ast }$.

Once Equation (2) is solved and $\theta \left(z\right)$ is determined, the normal stresses are evaluated via ${\sigma }_z\left(z,s\right)=-\frac{B\left(z\right)}{I_{\omega }}\omega \left(s\right)$, where $B\left(z\right)=E{I}_{\omega }{\theta }^{″}$ is the bimoment. Evaluation of the tangential stresses is more complicated. Indeed, they consist of an active and a reactive component, ${\tau }_{zs}={\tau }_{zs}^a+{\tau }_{zs}^r$, with the former being associated with the uniform flow of the shear strains and the latter being associated with the geometric constraint that enforces the uniformity. This circumstance differs from that occurring in open TWB, where tangential stresses are only of reactive nature. The active stress components are evaluated by the Bredt formula, ${\tau }_{zs}^a=\frac{M_t}{2{\mathsf{\Omega }}_{0}t}$. The reactive stress components are determined by integrating the indefinite equilibrium equation averaged on the thickness, ${\left({\tau }_{zs}^{r}t\right)}_{,s}+{\left({\sigma }_{z}t\right)}_{,z}=0$, where a comma denotes partial differentiation with respect the following variable. This supplies ${\tau }_{zs}^r\left(z,s\right)t\left(s\right)=-\frac{M_t^{\ast }\left(z\right)S_{\omega }\left(s\right)}{I_{\omega }}+q_0\left(z\right)$, where $S_{\omega }\left(s\right):={\oint }_{\Gamma }\omega \left(s\right)t\left(s\right)ds$ and $q_0\left(z\right)$ is an arbitrary function; the latter is determined by enforcing the static equivalence between ${\tau }_{zs}^r$ and $M_t^{\ast }$. Details on these calculations are given in Appendix A for rectangular box-girders.

2.2. GBT Theory

When the hypothesis of undeformability of the cross section is removed, the box-girder behaves as a plate assembly. According to the classic version of the GBT theory [23,37] (in which only the so-called conventional modes are taken into account), the displacement field is expressed as follows:

$$\begin{array}{c}u(z,s)=\sum _{k=1}^K{U}_k\left(s\right)\hspace{1em}{a}_k\left(z\right)={\mathit{U}}^T\mathbf{a}\hfill \\ v(z,s)=\sum _{k=1}^K{V}_k\left(s\right)\hspace{1em}{a}_k\left(z\right)={\mathit{V}}^T\mathbf{a}\hfill \\ w(z,s)=\sum _{k=1}^K{W}_k\left(s\right)\hspace{1em}{a}_k^{\prime }\left(z\right)={\mathit{W}}^T{\mathbf{a}}^{\prime }\hfill \end{array}$$

(3)

where $u,v,w$ are tangential (along the directrix), normal (along the thickness), and axial (along the axis) components, respectively. Here, $U_k\left(s\right),V_k\left(s\right),W_k\left(s\right)$ are assumed to be modal shapes and $a_k\left(z\right)$ are unknown amplitude functions. When the Kirchhoff hypothesis is invoked, the relevant strains $\mathit{\epsilon }={\left({\epsilon }_s^m,{\epsilon }_z^m,{\gamma }_{zs}^m,{\epsilon }_s^f,{\epsilon }_z^f,{\gamma }_{zs}^f\right)}^T$ are derived:

$$\begin{array}{c}\hfill {\epsilon }_s^m={\mathit{U}}^{\prime }{}^T\mathbf{a},\hspace{1em}{\epsilon }_z^m={\mathit{W}}^T{\mathbf{a}}^{″},\hspace{1em}{\gamma }_{zs}^m=\left({\mathit{U}}^T+{\mathit{W}}^{\prime }{}^T\right){\mathbf{a}}^{\prime }\\ \hfill {\epsilon }_s^f=-y{\mathit{V}}^{″}{}^T\mathbf{a},\hspace{1em}{\epsilon }_z^f=-y{\mathit{V}}^T{\mathbf{a}}^{″},\hspace{1em}{\gamma }_{zs}^f=-2y{\mathit{V}}^{\prime }{}^T{\mathbf{a}}^{\prime }\end{array}$$

(4)

where apexes $m,f$ denote membrane and flexural contributions, respectively; y is the distance from the middle plane; and a prime indicates differentiation with respect the independent variable. From linear Hooke’s law for isotropic material, the planar stresses $\mathit{\sigma }={\left({\sigma }_s^m,{\sigma }_z^m,{\tau }_{zs}^m,{\sigma }_s^f,{\sigma }_z^f,{\tau }_{zs}^f\right)}^T$ are found as follows:

$$\begin{array}{c}\hfill {\sigma }_s^m=E{\mathit{U}}^{\prime }{}^T\mathbf{a},\hspace{1em}{\sigma }_z^m=E{\mathit{W}}^T{\mathbf{a}}^{″},\hspace{1em}{\tau }_{zs}^m=G\left({\mathit{U}}^T+{\mathit{W}}^{\prime }{}^T\right){\mathbf{a}}^{\prime }\\ \hfill {\sigma }_s^f=-y\frac{E}{1-{\nu }^2}\left({\mathit{V}}^{″}{}^T\mathbf{a}+\nu {\mathit{V}}^T{\mathbf{a}}^{″}\right),\hspace{1em}{\sigma }_z^f=-y\frac{E}{1-{\nu }^2}\left(\nu {\mathit{V}}^{″}{}^T\mathbf{a}+{\mathit{V}}^T{\mathbf{a}}^{″}\right)\\ \hfill {\tau }_{zs}^f=-2yG{\mathit{V}}^{\prime }{}^T{\mathbf{a}}^{\prime }\end{array}$$

(5)

where $\nu$ is the Poisson ratio.

Finally, by making use of the Virtual Work Principle,

$$\underset{0}{\overset{\ell }{\int }}dz\underset{-\frac{t\left(s\right)}{2}}{\overset{\frac{t\left(s\right)}{2}}{\int }}dy\underset{\Gamma }{\oint }{\mathit{\sigma }}^T\delta \mathit{\epsilon }ds=\underset{0}{\overset{\ell }{\int }}dz\underset{\Gamma }{\oint }{\mathbf{f}}^T\delta \mathbf{u}ds$$

(6)

in which $\mathbf{f}={\left(f_s\left(z,s\right),f_y\left(z,s\right),f_z\left(z,s\right)\right)}^T$ are surface external forces and $\mathbf{u}=\left(u\right(z,s),v(z,s),$${w(z,s))}^T$, the following ordinary differential equations and relevant boundary conditions are obtained:

$$\begin{array}{cc}\hfill ({\mathit{C}}^e+{\mathit{C}}^f){\mathbf{a}}^{″″}+({\mathit{D}}^f+{\mathit{D}}^{fT}-{\mathit{D}}^s-{\mathit{D}}^t){\mathbf{a}}^{″}+({\mathit{B}}^f+{\mathit{B}}^d)\mathbf{a}& =\mathit{p}\hfill \\ \hfill {\left[\delta {\mathbf{a}}^{\prime }{}^T\left(\left({\mathit{C}}^f+{\mathit{C}}^e\right){\mathbf{a}}^{″}+{\mathit{D}}^f\mathbf{a}\right)\right]}_0^{\ell }& =\mathbf{0}\hfill \\ \hfill {\left[\delta {\mathbf{a}}^T\left(\left({\mathit{C}}^f+{\mathit{C}}^e\right){\mathbf{a}}^{‴}+\left({\mathit{D}}^f-{\mathit{D}}^s-{\mathit{D}}^t\right){\mathbf{a}}^{\prime }\right)\right]}_0\ell & =\mathbf{0}\hfill \end{array}$$

(7)

where

$$\begin{array}{c}{\mathit{C}}^e=\underset{\Gamma }{\oint }Et\mathit{W}{\mathit{W}}^{T}ds,\hspace{1em}{\mathit{C}}^f=\underset{\Gamma }{\oint }D\mathit{V}{\mathit{V}}^{T}ds\hfill \\ {\mathit{D}}^f=\underset{\Gamma }{\oint }\nu D\mathit{V}{\mathit{V}}^{″}{}^{T}ds,\hspace{1em}{\mathit{D}}^t=\underset{\Gamma }{\oint }\frac{G{t}^3}{3}{\mathit{V}}^{\prime }{\mathit{V}}^{\prime }{}^{T}ds,\hfill \\ {\mathit{D}}^s=\underset{\Gamma }{\oint }Gt{\left({\mathit{W}}^{\prime }{}^T+{\mathit{U}}^T\right)}^T\left({\mathit{W}}^{\prime }{}^T+{\mathit{U}}^T\right)ds\hfill \\ {\mathit{B}}^f=\underset{\Gamma }{\oint }D{\mathit{V}}^{″}{\mathit{V}}^{″}{}^{T}ds,\hspace{1em}{\mathit{B}}^d=\underset{\Gamma }{\oint }Et\mathit{U}{\mathit{U}}^{\prime }{}^{T}ds\hfill \end{array}$$

(8)

are known symmetric matrices, with $D\left(s\right):=\frac{E{t}^3\left(s\right)}{12\left(1-{\nu }^2\right)}$; moreover,

$$\mathit{p}=\underset{\Gamma }{\oint }{f}_s\mathit{U}ds+\underset{\Gamma }{\oint }{f}_y\mathit{V}ds-\underset{\Gamma }{\oint }\frac{\partial f_z}{\partial z}\mathit{W}ds$$

(9)

are generalized (modal) loads. Equation (7) is referred to as GBT equations. In them, matrices $\mathit{C},\mathit{D},\mathit{B}$ have been labeled as (f) flexural, (e) extensional, (s) shear, (t) torsional, and (d) in-plane dilatation, to remember their energy origin. In particular, ${\mathit{C}}^f$ accounts for flexure of the plates in the longitudinal plane, ${\mathit{D}}^f$ represents the transverse Poisson-induced flexure, and ${\mathit{B}}^f$ is thefor (direct) flexure in the transverse plane.

Once the GBT equations are solved and the amplitude functions $\mathbf{a}\left(z\right)$ are evaluated, the active stresses are computed from Equation (5). Reactive stresses, often ignored in the literature, must be determined by equilibrium arguments, as discussed below.

3. A Minimal GBT Model

A rectangular, $b\times h$, double-symmetric box-girder is considered (Figure 1). Its cross section is made of four thin elements: the bottom/upper flanges $i=1,3$, of thickness $t_1=t_3$, and the right/left webs $i=2,4$, of thickness $t_2=t_4$. A simple model is formulated, aimed at capturing the main mechanical behavior of the twist-distortional behavior. The task is accomplished by following the GBT approach, in which just two modes are considered: (i) the torsional mode, in which the cross section keeps its shape unaltered, and (ii) the distortional mode, in which the cross section is deformed. According to Equation (3), the local displacement field is expressed as follows:

$$\left(\begin{array}{c}u\left(z,s\right)\\ v\left(z,s\right)\\ w\left(z,s\right)\end{array}\right)=\theta \left(z\right)\left(\begin{array}{c}{U}_t\left(s\right)\\ V_t\left(s\right)\\ W_t\left(s\right)\end{array}\right)+\phi \left(z\right)\left(\begin{array}{c}{U}_d\left(s\right)\\ V_d\left(s\right)\\ W_d\left(s\right)\end{array}\right)$$

(10)

where the s-functions are the components of the two modes, to be properly chosen, and $\mathbf{a}={\left(\theta \left(z\right),\phi \left(z\right)\right)}^T$ are the unknown amplitude functions, measuring the magnitude of twist $\theta$ and distortion $\phi$.

To define the amplitudes, a general inextensional transformation of the box-girder is considered (by ignoring, at this stage, compatibility at corners), in which the x-axis rotates of an angle ${\psi }_x$ and the y-axis of an angle ${\psi }_y$, both positive counterclockwise (Figure 2a). By decomposing the transformation in its skew-symmetric (Figure 2b) and symmetric (Figure 2c) parts, the twist $\theta$ and the distortion $\phi$ are uniquely determined as follows:

$$\theta =\frac{1}{2}\left({\psi }_x+{\psi }_y\right),\hspace{1em}\phi =\frac{1}{2}\left({\psi }_x-{\psi }_y\right)$$

(11)

together with their inverse ${\psi }_x=\theta +\phi$, ${\psi }_y=\theta -\phi$.

To define the modes, a suitable in-plane field $U\left(s\right),V\left(s\right)$ is chosen, as illustrated ahead. From this, the warping component $W\left(s\right)$ is derived by exploiting the internal constraint condition ${\gamma }_{zs}\left(z,s\right)t\left(s\right)=\varphi \left(z\right)$, borrowed by the Bredt theory. Since, in a generic mode, ${\gamma }_{zs}=u_{,z}+w_{,s}=\left(U\left(s\right)+W^{\prime }\left(s\right)\right)a^{\prime }\left(z\right)$, the constraint entails $\left(U\left(s\right)+W^{\prime }\left(s\right)\right)t\left(s\right)=\mathrm{const}=:Q$. This equation can be integrated to furnish

$$W=-\underset{0}{\overset{s}{\int }}U\left(s\right)ds+Q\underset{0}{\overset{s}{\int }}\frac{ds}{t\left(s\right)}+C$$

(12)

where C is a further integration constant. The previous result is particularized to the two modes.

3.1. Torsional Mode

The in-plane displacement field is expressed as a rotation $\theta =1$ around the centroid, entailing (Figure 3a,b)

$$U_{t_i}=\left\{\begin{array}{cc}\frac{h}{2}\hfill & i=1,3\hfill \\ \frac{b}{2}\hfill & i=2,4\hfill \end{array}\right.,\hspace{1em}{V}_{t_i}=s_i\hspace{1em}i=1,\dots 4$$

(13)

where $s_i\hspace{1em}\left(i=1,\dots ,4\right)$ are local abscissas with origins at the midpoint of each element i. The associated warping is evaluated by Equation (12), applied to each element, which leads to the following:

$$W_{t_i}=\left\{\begin{array}{cc}\left(-\frac{h}{2}+\frac{Q_t}{t_i}\right)s_i+C_i\hfill & i=1,3\hfill \\ \left(-\frac{b}{2}+\frac{Q_t}{t_i}\right)s_i+C_i\hfill & i=2,4\hfill \end{array}\right.$$

(14)

where $C_i$ are integration constants. By enforcing continuity at the four corners and the condition ${\oint }_{\Gamma }W\left(s\right)t\left(s\right)ds=0$ (for equilibrium along z), the five unknown are determined as $C_i=0$ and $Q_t=\frac{bh{t}_1{t}_2}{b{t}_2+h{t}_1}$. Therefore,

$$W_{t_i}=\left\{\begin{array}{cc}\frac{2w_t}{b}{s}_i\hfill & i=1,3\hfill \\ -\frac{2w_t}{h}{s}_i\hfill & i=2,4\hfill \end{array}\right.$$

(15)

where $w_t:=\frac{bh}{4}\left(\frac{b{t}_2-h{t}_1}{b{t}_2+h{t}_1}\right)$ is the modulus of warping at the corners (Figure 3c). This expression coincides with that in Equation (A1) of Appendix A, obtained by the Vlasov theory; therefore, $W_{t_i}\equiv \omega \left(s\right)$. It should be noticed that the torsional warping vanishes when $\frac{b}{h}=\frac{t_1}{t_2}.$

3.2. Distortional Mode

The in-plane displacement field is found as that experienced by a planar frame, undergoing displacements at joints, namely (i) prescribed translation of joints, caused by a distortion $\phi =1$, with no rotations, (Figure 4a); (ii) unknown rotations of joints of equal amplitude $\alpha$ with no translations (Figure 4b). After having determined the moments at any of the the four joints, equilibrium is enforced, from which the unknown is evaluated as follows:

$$\alpha :=-\frac{b{t}_2^3-h{t}_1^3}{b{t}_2^3+h{t}_1^3}$$

(16)

referred to below as the distortional joint rotation (which can be positive, nil, or negative according to the geometric parameters. After that, the transverse displacements are computed by integrating the elastic line equation under now-known translations and rotations at the ends (see Appendix B for details). The tangential displacements are instead immediately evaluated from the kinematics illustrated in Figure 4a. By summarizing (Figure 5a,b),

$$U_{d_i}=\left\{\begin{array}{cc}-\frac{h}{2}\hfill & i=1,3\hfill \\ \frac{b}{2}\hfill & i=2,4\hfill \end{array}\right.,\hspace{1em}{V}_{d_i}=\left\{\begin{array}{cc}\frac{1}{2}{s}_i\left(3+\frac{4s_i^2\left(\alpha -1\right)}{b^2}-\alpha \right)\hfill & i=1,3\hfill \\ \frac{1}{2}{s}_i\left(-3+\frac{4s_i^2\left(\alpha +1\right)}{h^2}-\alpha \right)\hfill & i=2,4\hfill \end{array}\right.$$

(17)

to obtain the warping component, the procedure already described for the torsional mode is followed, leading to the following:

$$W_{d_i}=\left\{\begin{array}{cc}\left(\frac{h}{2}+\frac{Q_d}{t_i}\right)s_i+C_i\hfill & i=1,3\hfill \\ \left(-\frac{b}{2}+\frac{Q_d}{t_i}\right)s_i+C_i\hfill & i=2,4\hfill \end{array}\right.$$

(18)

Due to the skew-symmetric nature of $U_{d_i},$ continuity at joints and the zero-average condition require $C_i=Q_d=0$. Therefore,

$$W_{d_i}=\left\{\begin{array}{cc}\frac{2w_d}{b}{s}_i\hfill & i=1,3\hfill \\ -\frac{2w_d}{h}{s}_i\hfill & i=2,4\hfill \end{array}\right.$$

(19)

where $w_d:=\frac{bh}{4}$ is the modulus of the warping at joints (Figure 5c).

It is worth noticing that, differently from torsion, the distortional warping never vanishes. The two warping components are proportional (Figure 3c and Figure 5c), i.e., $W_t\left(s\right)=\beta W_d\left(s\right)$, with

$$\beta :=\frac{b{t}_2-h{t}_1}{b{t}_2+h{t}_1}$$

(20)

defined as the warping ratio. The plot of $\beta$ vs. the aspect ratio $\frac{b}{h}$, for different thickness ratios $\frac{t_1}{t_2}$, is represented in Figure 6. In the range examined, $\beta$ spans the interval $\left(-1,1\right),$ so that the two warpings (for unitary twist and distortion) are of the same order of magnitude, except close to a critical combination of the parameters, for which torsional warping is zero.

3.3. GBT Equations

By using the torsional mode (Equations (13) and (15)) and the distortional mode (Equations (17) and (19)), the $2\times 2$ matrices (Equation (8) appearing in the GBT equations (Equation (7)) are evaluated, for which the components are as follows:

$$\begin{array}{c}{C}_{\alpha \beta }^e=E\underset{\Gamma }{\oint }t{W}_{\alpha }{W}_{\beta }ds,\hspace{1em}{C}_{\alpha \beta }^f=\underset{\Gamma }{\oint }D{V}_{\alpha }{V}_{\beta }ds\hfill \\ D_{\alpha \beta }^f=\nu \underset{\Gamma }{\oint }{V}_{\alpha }{V}_{\beta }^{″}ds,\hspace{1em}{D}_{\alpha \beta }^t=\frac{G}{3}\underset{\Gamma }{\oint }{t}^3{V}_{\alpha }^{\prime }{V}_{\beta }^{\prime }ds\hfill \\ D_{\alpha \beta }^s=G\underset{\Gamma }{\oint }t(W_{\alpha }^{\prime }+U_{\alpha })(W_{\beta }^{\prime }+U_{\beta })ds\hfill \\ B_{\alpha \beta }^f=\frac{E}{12}\underset{\Gamma }{\oint }{t}^3{V}_{\alpha }^{″}{V}_{\beta }^{″}ds,\hspace{1em}{B}_{\alpha \beta }^d=E\underset{\Gamma }{\oint }t{U}_{\alpha }{U}_{\beta }^{\prime }ds\hfill \\ \alpha ,\beta =t,d\hfill \end{array}$$

(21)

together with $\mathbf{a}={\left(\theta ,\phi \right)}^T$, $\mathit{p}={\left(c_t,c_d\right)}^T$, in which $c_{\alpha }$ is the couples per unit length. A careful check of these matrices reveals the following:

• since $U_t\left(s\right)$ and $U_d\left(s\right)$ are step-wise constant, $U_t^{\prime }=U_d^{\prime }=0$, and therefore ${\mathit{B}}^d=\mathbf{0}$ (i.e., no in-plane dilatation occurs);
• since $V_t\left(s\right)$ is step-wise linear and $V_d\left(s\right)$ quadratic, then $B_{11}^f=B_{12}^f=B_{21}^f=0$, so that only $B_{22}^f\ne 0$ in matrix ${\mathit{B}}^f$;
• since $W_t^{\prime }+U_t=Q_t\ne 0$ (torsional shear different from zero) and $W_d^{\prime }+U_d=Q_d=0$ (distortional shear zero), then $D_{12}^s=D_{21}^s=D_{22}^s=0$, so that only $D_{11}^s\ne 0$ in matrix ${\mathit{D}}^s$; the implications on stresses of $Q_d=0$ will be discussed ahead.

Moreover, aimed at obtaining the simplest model, the Poisson ratio is taken as zero, so that ${\mathit{D}}^f=\mathbf{0}$. Concerning ${\mathit{D}}^t$, it accounts for linear distribution of the tangential stresses in the thickness. When TWBs are considered, these terms are usually negligible with respect to the Bredt stresses; however, in the case of reinforced concrete box-girder bridges, where thicknesses are not so small, they can give significant contributions, as it will be shown by numerical simulations. Therefore, matrix ${\mathit{D}}^t$ is retained in the analysis.

As a result, the field GBT equations reduce to the following:

$$\left(\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right)\left(\begin{array}{c}{\theta }^{″″}\\ {\phi }^{″″}\end{array}\right)+\left(\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right)\left(\begin{array}{c}{\theta }^{″}\\ {\phi }^{″}\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& B_{22}^f\end{array}\right)\left(\begin{array}{c}\theta \\ \phi \end{array}\right)=\left(\begin{array}{c}{c}_t\\ c_d\end{array}\right)$$

(22)

where $C_{ij}:=C_{ij}^e+C_{ij}^f$ and $D_{ij}:=D_{ij}^s+D_{ij}^t$, together with the following boundary conditions:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{cc}\hfill {\left[{\left(\begin{array}{cc}\delta {\theta }^{\prime }& \delta {\phi }^{\prime }\end{array}\right)}^T\left(\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right)\left(\begin{array}{c}{\theta }^{″}\\ {\phi }^{″}\end{array}\right)\right]}_0^{\ell }& =\left(\begin{array}{c}0\\ 0\end{array}\right)\hfill \\ \hfill {\left[{\left(\begin{array}{cc}\delta \theta & \delta \phi \end{array}\right)}^T\left(\left(\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right)\left(\begin{array}{c}{\theta }^{‴}\\ {\phi }^{‴}\end{array}\right)-\left(\begin{array}{cc}{D}_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right)\left(\begin{array}{c}{\theta }^{\prime }\\ {\phi }^{\prime }\end{array}\right)\right)\right]}_0^{\ell }& =\left(\begin{array}{c}0\\ 0\end{array}\right)\hfill \end{array}\end{array}\end{array}$$

(23)

The explicit expressions of the matrices appearing in Equation (22) are reported in Appendix C. It is worth noticing that, when the distortion $\phi$ is ignored and the flexure of plates is neglected with respect to the membrane strains (i.e., $C_{11}^f\simeq 0,\hspace{1em}{D}_{11}^t\simeq 0$), the GBT equations reduce to Equation (2) of the Vlasov theory, with $C_{11}^e=E{I}_{\omega }$, $D_{11}^s=GJ$. When torsion is ignored (together with $C_{22}^f\simeq 0,\hspace{1em}{D}_{22}^t\simeq 0$) and distortion is artificially uncoupled, the GBT equations reduce to the classic equation of the beam on elastic soil, on which the BEF analogy is founded.

3.4. Stresses

Once the GBT equations are solved and the amplitudes $\theta \left(z\right),\phi \left(z\right)$ are evaluated, the active stresses can be computed by coming back to Equation (5). Reactive stresses, instead, must be determined from equilibrium. By separating the effects of the two modes, the following results are drawn (with the Poisson ratio taken zero). Moreover, from now on, $\tau :={\tau }_{zs}$ is used to simplify the notation.

In the torsional mode, the active stresses are as follows:

$$\begin{array}{c}{\sigma }_s^m(z,s)=0\hfill \\ {\sigma }_z^m(z,s)=\left\{\begin{array}{cc}\frac{E}{2}h\beta s_i{\theta }^{″}\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ -\frac{E}{2}b\beta s_i{\theta }^{″}\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \\ {\tau }^m(z,s)=\left\{\begin{array}{cc}G\frac{bh{t}_2}{h{t}_1+b{t}_2}{\theta }^{\prime }\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ G\frac{bh{t}_1}{h{t}_1+b{t}_2}{\theta }^{\prime }\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \\ {\sigma }_s^f(z,s)=0\hfill \\ {\sigma }_z^f(z,s)=-E{y}_i{s}_i{\theta }^{″}\left(z\right)\hspace{1em}i=1,\dots 4\hfill \\ {\tau }^f(z,s)=-2G{y}_i{\theta }^{\prime }\left(z\right)\hspace{1em}i=1,\dots 4\hfill \end{array}$$

(24)

and the reactive stresses (see the Appendix A) are as follows:

$${\tau }^{mr}\left(z,s\right)=\left\{\begin{array}{cc}-E{\theta }^{‴}\left(z\right)h\frac{\left(h{t}_1-b{t}_2\right)\left(3h\left(b^2-4s_i^2\right)t_1+b\left(b^2+2\left(h^2-6s_i^2\right)\right)t_2\right)}{48{\left(h{t}_1+b{t}_2\right)}^2}\hfill & i=1,3\hfill \\ E{\theta }^{‴}\left(z\right)b\frac{\left(h{t}_1-b{t}_2\right)\left(h\left(2b^2+h^2-12s_i^2\right)t_1+3b\left(h^2-4s_i^2\right)t_2\right)}{48{\left(h{t}_1+b{t}_2\right)}^2}\hfill & i=2,4\hfill \end{array}\right.$$

(25)

The normal stresses ${\sigma }_z^m(z,s)$ coincide with those predicted by the Vlasov theory (see Equation (A3) in Appendix A), while the tangential stresses ${\tau }^m(z,s)$ coincide with those of the Bredt theory (see Equation (A4) in Appendix A). The reactive tangential stresses ${\tau }^{mr}$ cannot be captured by the constitutive law of the GBT model (which is based on the same Vlasov kinematic constraint on the flow).

The relevant stress-fields are qualitatively depicted in Figure 7 at a generic cross-section; there, the index max denotes evaluation at the top/bottom face of the plates (i.e., at $y_i=\pm \frac{t_i}{2}$). The diagrams repeat themselves at any abscissas z, being altered by a factor proportional to ${\theta }^{\prime }\left(z\right),{\theta }^{″}\left(z\right),{\theta }^{‴}\left(z\right)$, as stated by Equations (24) and (25). The membrane normal stresses ${\sigma }_z^m$ (Figure 7a) are caused by nonuniform warping; they are proportional to $W_t\left(s\right)$ (Equation (15)) and vary on z as ${\theta }^{″}\left(z\right)$. The flexural normal stresses ${\sigma }_z^f$ (maximum value in Figure 7b) are caused by the flexure of plates in the longitudinal direction; they are proportional to $V_t\left(s\right)$ (Equation (13)) and depend on z via ${\theta }^{″}\left(z\right)$. The membrane active tangential stresses ${\tau }^m$ (Figure 7c) are consistent with the Bredt behavior of the box; they are step-wise constant on the cross section and vary with ${\theta }^{\prime }\left(z\right)$. The reactive component ${\tau }^{mr}$ (Figure 7d) must be added to them; they are step-wise parabolic on the cross section and vary with ${\theta }^{‴}\left(z\right)$. The flexural tangential stresses ${\tau }^f$ (maximum value in Figure 7e) are related to the torsion of the plates; being proportional to $V_t^{\prime }\left(s\right)$, they are step-wise constant on the cross section and vary with ${\theta }^{\prime }\left(z\right)$.

In the distortional mode, the active stresses are as follows:

$$\begin{array}{c}{\sigma }_s^m(s,z)=0\hfill \\ {\sigma }_z^m(s,z)=\left\{\begin{array}{cc}\frac{E}{2}h{s}_i{\phi }^{″}\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ -\frac{E}{2}b{s}_i{\phi }^{″}\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \\ {\tau }^m(s,z)=0\hfill \\ {\sigma }_s^f(s,z)=\left\{\begin{array}{cc}-E{y}_i\frac{12\left(\alpha -1\right)}{b^2}{s}_i\phi \left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ -E{y}_i\frac{12\left(\alpha +1\right)}{h^2}{s}_i\phi \left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \\ {\sigma }_z^f(s,z)=\left\{\begin{array}{cc}-\frac{E}{2}{y}_i{s}_i\left(3+4\left(\alpha -1\right)\frac{s_i^2}{b^2}-\alpha \right){\phi }^{″}\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ -\frac{E}{2}{y}_i{s}_i\left(-3+4\left(\alpha +1\right)\frac{s_i^2}{h^2}-\alpha \right){\phi }^{″}\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \\ {\tau }^f(s,z)=\left\{\begin{array}{cc}G\frac{y_i}{b^2}\left(b^2\left(\alpha -3\right)-12\left(\alpha -1\right)s_i^2\right){\phi }^{\prime }\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ G\frac{y_i}{h^2}\left(h^2\left(\alpha +3\right)-12s_i^2\left(\alpha +1\right)\right){\phi }^{\prime }\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.\hfill \end{array}$$

(26)

with $\alpha$ defined in Equation (16). For the torsional mode as well, the membrane reactive tangential stresses must be determined from equilibrium, leading to the following (see the Appendix D):

$${\tau }^{mr}\left(z,s\right)=\left\{\begin{array}{cc}\left(\frac{bh}{48}\frac{2b{t}_1+h{t}_2}{t_1}-\frac{1}{4}h{s}_i^2\right)E{\phi }^{‴}\left(z\right)\hfill & i=1,\hspace{1em}3\hfill \\ \left(-\frac{bh}{48}\frac{\left(b{t}_1+2h{t}_2\right)}{t_2}+\frac{1}{4}b{s}_i^2\right)E{\phi }^{‴}\left(z\right)\hfill & i=2,\hspace{1em}4\hfill \end{array}\right.$$

(27)

The relevant stress-fields (with max denoting evaluation at $y_i=\pm \frac{t_i}{2}$) are qualitatively depicted in Figure 8. The diagrams repeat themselves at any abscissas z, being altered by a factor proportional to $\phi \left(z\right),{\phi }^{\prime }\left(z\right),{\phi }^{″}\left(z\right)$, as stated by Equations (26) and (27). The membrane normal stresses ${\sigma }_z^m$ (Figure 8a) are caused by nonuniform warping; they are proportional to $W_d\left(s\right)$ (Equation (19)) and vary on z as ${\phi }^{″}\left(z\right)$. The flexural normal stresses ${\sigma }_z^f$ (maximum value in Figure 8b) are due to the longitudinal flexure of plates; they, being proportional to $V_d\left(s\right)$ (Equation (17)), are cubic on the cross section and depend on z via ${\phi }^{″}\left(z\right)$. The transverse flexural normal stresses ${\sigma }_s^f$ (maximum value in Figure 8c) are generated by the flexure of the plates in the transverse direction; being proportional to $V_d^{″}\left(s\right)$, they are step-wise linear on the cross section and are modulated by $\phi \left(z\right)$. The membrane reactive tangential stresses ${\tau }^{mr}$ (Figure 8d) are step-wise parabolic and proportional to ${\phi }^{‴}\left(z\right)$. The flexural tangential stresses ${\tau }^f$ (maximum value in Figure 8e) are related to the torsion of the plates; they depend on $V_d^{\prime }\left(s\right)$, are step-wise parabolic on the cross section, and are modulated by ${\phi }^{\prime }\left(z\right)$.

It is worth noticing that, according to the model presented here, the distortional tangential stresses are purely reactive, since no shear strains occur in the distortional mode; in contrast, longitudinal normal stresses are of the active type, directly related to distortional warping. The model, therefore, is different from that adopted in Reference [8], where due to the assumed step-wise constant shear with no warping, the tangential stress is active while the normal stress is reactive. As an advantage of the method followed here, the continuity of the flow of tangential stresses is assured at the joints, a condition that cannot be satisfied in the alternative approach.

4. Algorithmic Aspects

Aimed at further simplifying the reduced GBT equations (Equation (22)), an analysis of the relative magnitude of different terms is carried out. Moreover, qualitative information on the relative importance between distortion and twist are sought. Then, some methods of solution are discussed.

4.1. Order of Magnitude of the Coefficients and Unknowns

The possibility to neglect some terms in GBT equations depending on the geometric quantities involved is detected.

4.1.1. Flexural vs. Extensional Higher-Order Derivatives

First, the coefficients of the higher-order derivatives, $C_{ij}=C_{ij}^e+C_{ij}^f$, are considered. In them, it is expected that the $C_{ij}^f$s, which account for longitudinal flexure of the thin plates, are much smaller than the $C_{ij}^e$s, which account for extension of the longitudinal fibers induced by nonuniform warping (at least in the distortional mode, where warping never disappears). To investigate this aspect, the ratios $C_{ij}^e/C_{ij}^f$ are plotted in Figure 9 vs. the $\frac{h}{t}$ ratio for different $\frac{b}{h}$ ratios and for a fixed $\frac{t_2}{t_1}$. It is seen that, for realistic values of $\frac{h}{t}$, e.g., ranging between 5 and 10, the flexural coefficients are, as expected, a small fraction (few percents) of the extensional ones. It is concluded that the $C_{ij}^f$ are generally negligible, especially for large $\frac{h}{t}$. An exception, of course, occurs when $C_{11}^e=C_{12}^e=0$.

4.1.2. Extensional vs. Shear Torsional Effects

A second, and more important, topic concerns coupling existing between the GBT equations. They are a set of two fourth-order differential equations, where coupling is mainly due to higher-order terms (the larger ones being produced by warping, as seen before). It is interesting to investigate if, and under which conditions, such a coupling can be neglected. To this end, the box-girder is assumed to be loaded by sinusoidal loads $\mathit{p}={\left(\widehat{c_t},{\widehat{c}}_d\right)}^{T}sin\left(\frac{n\pi z}{\ell }\right)$, which trigger sinusoidal deflections $\mathbf{a}={\left(\widehat{\theta },\widehat{\phi }\right)}^{T}sin\left(\frac{n\pi z}{\ell }\right)$. The amplitudes must satisfy the following algebraic equations:

$$\left[\left(\begin{array}{cc}{C}_{11}^e& C_{12}^e\\ C_{21}^e& C_{22}^e\end{array}\right){\left(\frac{n\pi }{\ell }\right)}^4+\left(\begin{array}{cc}{D}_{11}^s& 0\\ 0& 0\end{array}\right){\left(\frac{n\pi }{\ell }\right)}^2+\left(\begin{array}{cc}0& 0\\ 0& B_{22}^f\end{array}\right)\right]\left(\begin{array}{c}\widehat{\theta }\\ \widehat{\phi }\end{array}\right)=\widehat{c}\left(\begin{array}{c}1\\ 1\end{array}\right)$$

(28)

where $\widehat{c_t}={\widehat{c}}_d=\widehat{c}$ has been taken and $C_{ij}^f$ is neglected together with $D_{ij}^t$ (small thickness). The task is understanding if the $C_{ij}^e$s bring a significant contribution to the equilibrium when, e.g., they are compared to $D_{11}^s=GJ$ or $B_{22}^f$. Referring to a rectangular box-girder with uniform thickness $t_1=t_2=:t$, the nondimensional ratios $r_{ij}^s:=\frac{C_{ij}^e}{D_{11}^s}{\left(\frac{n\pi }{\ell }\right)}^2$ and $r_{ij}^f:=\frac{C_{ij}^e}{B_{22}^f}{\left(\frac{n\pi }{\ell }\right)}^4$ among the coefficients are plotted in Figure 10 vs. the squatness ratio $\frac{nh}{\ell }$ for different $\frac{b}{h}$. It is seen that the $r_{ij}^s$ are less than 0.1 if $\frac{nh}{\ell }<0.2$ (e.g., $n<4$ if $\ell =20 h$). Moreover, $r_{ij}^s$ is much smaller as $\frac{b}{h}$ decreases, tending towards a square shape. Concerning $r_{ij}^f$, similar consideration hold.

The investigation therefore leads to the following conclusions:

• Warping, as expected, plays a minor but not negligible role in torsion of closed TWB. In particular, it is fundamental in describing boundary layers close the constraints (or lumped forces), which, in a Fourier perspective, call for higher-order harmonics (n large), which make the fourth-order derivatives comparable with second- or zero-order derivatives.
• Due to the small but not negligible coupling terms due to warping, the torsion–distortion mechanical problem cannot, in principle, be uncoupled, as already observed, e.g., in Reference [8]. A measure of the error made in splitting the problem will be discussed with reference to the numerical results.

4.1.3. Distortional vs. Twist Amplitude

If, in a first attempt, the warping effects are completely neglected (i.e., $C_{ij}^e=0$ is taken) together with the torsional curvatures (i.e., $D_{ij}^t=0$), the two Equation (28) uncouple. Their solution gives a rough estimation of the relative importance of the distortion on twist, namely

$$\frac{\widehat{\phi }}{\widehat{\theta }}={\left(\frac{\pi }{\ell }\right)}^2\frac{D_{11}^s}{B_{22}^f}$$

(29)

This ratio is plotted in Figure 11 vs. $\frac{t}{h}$ for squared (Figure 11a) and rectangular (Figure 11b) boxed cross sections having uniform thickness $t_1=t_2=:t$ and for different slenderness ratios $\frac{\ell }{nh}$. It appears that distortion is of the same order as twist or even larger:

• for a fixed slenderness ratio $\frac{\ell }{nh}$, distortion is larger for smaller thicknesses;
• for a fixed thickness ratio $\frac{t}{h}$, distortion is larger for shorter lengths.

The results corroborate and quantify the common idea that the girder behaves as a beam (i.e., with no distortion) when it is long and its cross section is thick.

4.2. Solution Methods

The exact integration of the complete GBT equations, although not difficult in principle, is quite laborious. As a matter of fact, when boundary layers exist, even exact closed-form solution are not well-conditioned and numerical problems arise. To confine the analysis to a minimum level, two approaches are proposed here: (a) an exact Fourier analysis, holding for simply supported box-girder with free warping at the ends; (b) an exact closed-form solution for simplified equations, heuristically obtained by neglecting coupling.

4.2.1. Fourier Analysis

The modal loads $c_t,c_d$ are expanded in Fourier sinus series as $c_{\alpha }={\sum }_{n=1}^N{c}_{{\alpha }_n}sin\left(\frac{n\pi }{\ell }z\right)$ with $\alpha =t,d$. Accordingly, generalized displacements are expanded as well as

$$\left(\begin{array}{c}\theta \\ \phi \end{array}\right)=\sum _{n=1}^N\left(\begin{array}{c}{\theta }_n\\ {\phi }_n\end{array}\right)sin\left(\frac{n\pi }{\ell }z\right)$$

(30)

with ${\theta }_n,{\phi }_n$ unknowns. By substituting the series in the GBT equations (Equation (22)) and by separating the harmonics, we obtain the following (apex omitted on $B_{22}$):

$$\left(\begin{array}{cc}\frac{D_{11}{\pi }^2{n}^2}{{\ell }^2}+\frac{C_{11}{\pi }^4{n}^4}{{\ell }^4}& \frac{D_{12}{\pi }^2{n}^2}{{\ell }^2}+\frac{C_{12}{\pi }^4{n}^4}{{\ell }^4}\\ \frac{D_{21}{\pi }^2{n}^2}{{\ell }^2}+\frac{C_{21}{\pi }^4{n}^4}{{\ell }^4}& B_{22}+\frac{D_{22}{\pi }^2{n}^2}{{\ell }^2}+\frac{C_{22}{\pi }^4{n}^4}{{\ell }^4}\end{array}\right)\left(\begin{array}{c}{\theta }_n\\ {\phi }_n\end{array}\right)=\left(\begin{array}{c}{c}_{t_n}\\ c_{d_n}\end{array}\right)$$

(31)

From these equations, the unknowns are evaluated.

4.2.2. A Simplified Approach: The Uncoupled Equations

A heuristic approach is attempted here, consistent with the literature, in which all out-of-diagonal terms are (arbitrarily) neglected. In this case, the differential equations and boundary conditions read as follows:

$$\begin{array}{cc}\hfill E{I}_{\omega }{\theta }^{″″}-GJ{\theta }^{″}& =p_t\hfill \\ \hfill \theta \left(0\right) ={\theta }^{″}\left(0\right) =\theta \left(\ell \right) ={\theta }^{″}\left(\ell \right)& =0\hfill \end{array}$$

(32)

together with

$$\begin{array}{cc}\hfill C_{22}{\phi }^{″″}+B_{22}\phi & =p_d\hfill \\ \hfill \phi \left(0\right) ={\phi }^{″}\left(0\right) =\phi \left(\ell \right) ={\phi }^{″}\left(\ell \right)& =0\hfill \end{array}$$

(33)

The first equation coincides with Equation (2) of the Vlasov theory, while the second one is the classic equation of a beam on Winkler soil. It should be noticed that, according to this simplified approach, distortional stresses simply add themselves to Vlasov stresses without any feedback effect.

Integration usually calls for dividing the $\left(0,\ell \right)$ interval in subintervals in which the loads are continuous; in the generic subinterval i, the general solution reads as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\theta \left(z_i\right) =c_{1_i}{e}^{\lambda \left(z_i-l_i\right)}+c_{2_i}{e}^{-\lambda z_i}+c_{3_i}{z}_i+c_{4_i}+{\overline{\theta }}_i\left(z_i\right)\hfill \\ \phi \left(z_i\right) =e^{-\mu z}\left(c_{1_i}sin\left(\mu z\right) + c_{2_i}cos\left(\mu z\right)\right) + e^{\mu z}\left(c_{3_i}sin\left(\mu z\right) + c_{4_i}cos\left(\mu z\right)\right) + {\overline{\phi }}_i\left(z_i\right)\hfill \end{array}\end{array}$$

(34)

where $l_i$ is the lengths of the subintervals, $z_i$ is the local abscissas spanning them, with the origin at the left end; and $c_{k_i}$ is an arbitrary constant. Moreover, ${\overline{\theta }}_i\left(z_i\right),\hspace{1em}{\overline{\phi }}_i\left(z_i\right)$ are particular solutions. Finally,

$$\lambda :=\sqrt{\frac{GJ}{E{I}_{\omega }}},\hspace{1em}\mu :=\sqrt[4]{\frac{B_{22}}{4\hspace{1em}{C}_{22}}}$$

(35)

are wave numbers.

The arbitrary constants $c_j^{\left(i\right)}$ appearing in each problem are determined by enforcing two boundary conditions at each $z=0,\ell$ end and four conditions at each of the internal boundaries, as supplied by Equation (23). If no lumped couples are present, these read as follows:

• for the twist problem, (i) continuity of $\theta$, (ii) continuity of ${\theta }^{\prime }$, (iii) equilibrium of forces dual of $\delta \theta$ (entailing continuity of the bimoment $B=E{I}_{\omega }{\theta }^{″}$), and (iv) equilibrium of forces dual of $\delta {\theta }^{\prime }$ (entailing continuity of the total torsional moment ${\overline{M}}_t=GJ{\theta }^{\prime }-E{I}_{\omega }{\theta }^{‴}$);
• for the distortional problem, (i) continuity of $\phi$, (ii) continuity of ${\phi }^{\prime }$, (iii) equilibrium of forces dual of $\delta \phi$ (entailing continuity of ${\phi }^{″}$), and (iv) equilibrium of forces dual of $\delta {\phi }^{\prime }$ (entailing continuity of ${\phi }^{‴}$).

By numerically solving two sets of linear algebraic problems, the arbitrary constants are evaluated and the solution (Equation (34)) is analytically obtained in step-wise form.

5. Numerical Results

A case study is considered, consisting of a reinforced concrete, rectangular box-girder bridge, taken from the literature [3]. The box-girder has the following geometric characteristics: length $\ell =30 \mathrm{m}$, width $b=6 \mathrm{m}$, height $h=1.5 \mathrm{m}$, flange thickness $t_1=0.25 \mathrm{m}$, web thickness $t_2=0.35 \mathrm{m}$; elastic moduli $E=\mathrm{35,654}\times {10}^3 \mathrm{kN}/{\mathrm{m}}^2$,and $G=E/2$. The GBT constants consequently assume the values in

Table 1. Coefficients of the Generalized Beam Theory (GBT) equations.

(i,j)	(1,1)	(1,2)	(2,2)
$C_{ij}^e$	$1.18\times {10}^8$	$1.70\times {10}^8$	$2.44\times {10}^8$
$C_{ij}^f$	$1.74\times {10}^6$	$2.21\times {10}^6$	$3.29\times {10}^6$
$D_{ij}^s$	$1.02\times {10}^8$	0	0
$D_{ij}^t$	$1.88\times {10}^6$	$3.50\times {10}^5$	$2.63\times {10}^6$
$B_{ij}^f$	0	0	$6.80\times {10}^5$

. Moreover, the joint distortional rotation is $\alpha =-0.83$ and the warping ratio $\beta =0.70$.

The girder is assumed to be simply supported at the end, where it is free to warp. It is loaded at the two webs by equal and opposite forces $P=\pm 750 \mathrm{kN}$, uniformly distributed on a length $\mathsf{\Delta }=\frac{\ell }{4}$, as $p_0:=\pm \frac{P}{\mathsf{\Delta }}=100 \mathrm{kN}/\mathrm{m}$, and centered at the abscissa: $\xi =\frac{\ell }{4}$ (Figure 12). The relevant modal loads $c_t\left(z\right),c_d\left(z\right)$ are evaluated by Equation (9), as the virtual work spent at the abscissa z by the external forces in the displacements caused by $\theta =1$ and $\phi =1$, respectively. Since, at the application points of the loads (i.e., at the corners $C,D$ of the box, see Figure 1), it is $V=\pm \frac{b}{2}$ in both modes, it follows that

$$c_{\alpha }\left(x\right)=\left\{\begin{array}{cc}{p}_{0}b\hfill & z\in \left(\xi -\frac{\mathsf{\Delta }}{2},\xi +\frac{\mathsf{\Delta }}{2}\right)\hfill \\ 0\hfill & z\notin \left(\xi -\frac{\mathsf{\Delta }}{2},\xi +\frac{\mathsf{\Delta }}{2}\right)\hfill \end{array}\right.\hspace{1em}\alpha =t,d$$

(36)

The elastic problem was tackled in different ways: (a) by solving the GBT Equation (22) (with no terms neglected) via Fourier series; (b) by exactly integrating the differential Equations (32) and (33) after having neglected the coupling terms; and (c) via a FEM analysis. The latter was carried out by using a commercial software by implementing a model made of 40 × 60 = 2400 shell elements (along the directrix and the axis, respectively) for an overall number of 2440 degrees of freedom. The results were then compared.

5.1. Deflection Analysis

First, Fourier analysis was carried out for the GBT equations. By expanding the modal loads in sinus Fourier series, the coefficients $c_{t_n}=c_{d_n}=:c_n$ were found, as reported in Figure 13a, together with the truncated series with $N=20,50,100$ terms. The Fourier spectra for deflections, as furnished by Equation (31), together with the reconstructed series (for $N=50$, value used from now on) are represented in Figure 13b,c. As a general comment, distortion of the box girder is remarkable, larger than twist, and it propagates to a large distance from the loaded region.

The simplified method was successively applied by uncoupling the GBT equations. The integration domain was broken in three subintervals, i.e., ${\mathcal{I}}_1:=\left(0,\xi -\frac{\mathsf{\Delta }}{2}\right), {\mathcal{I}}_2:=\left(\xi -\frac{\mathsf{\Delta }}{2},\xi +\frac{\mathsf{\Delta }}{2}\right), {\mathcal{I}}_3:=\left(\xi +\frac{\mathsf{\Delta }}{2},\ell \right)$, and the general solution (Equation (34)) was written in each of them. There, ${\overline{\theta }}_i\left(z_i\right)=\overline{{\phi }_i}\left(z_i\right)=0 \left(i=1,3\right)$ and ${\overline{\theta }}_2\left(z_2\right)=-\frac{b{p}_0}{2GJ}{z}_2^2, {\overline{\phi }}_2\left(z_2\right)=\frac{b{p}_0}{B_{22}}$ are the particular solutions. The twelve arbitrary constants $c_{k_i}$ appearing in each problem were determined by enforcing two boundary conditions at each $z=0,\ell$ end and four conditions at each of the internal boundaries, as discussed before. The responses $\theta \left(z\right),\phi \left(z\right)$ thus obtained, Equation (34), are plotted in Figure 14 and compared with those of the complete Fourier analysis. Here, the exact numerical results, as provided by the FEM analysis are also reported. It appears that FEM and Fourier analysis are in excellent agreement, in spite of the stronger simplifications introduced in the analytical model. The simplified analysis also supplies satisfactory results (maximum error of about +11% on twist and about +10% on distortion), especially far from the loads. Differences with the exact results put into light the contribution of coupling, which should not be ignored, if more accurate results are desired. It is concluded that the error provided by the simplified approach is reasonable and that, at least in this case study, this method is conservative since it magnifies the deflection (in agreement with the conclusions of Reference [6]).

The first, second, and third derivatives of both $\theta \left(z\right)$ and $\phi \left(z\right)$ amplitudes, needed for stresses evaluation according to Equations (24)–(27), were then computed. In Figure 15, the Fourier (coupled) solution is compared with that of the simplified (uncoupled) method (comparison with FEM being not straightforward). It appears that, although small spurious oscillations persist with $N=50$, the Fourier solution is quite smooth. Concerning the third derivatives ${\theta }^{‴},{\phi }^{‴}$, two peaks are noticed at the end-points of the loaded interval, denoting the occurrence of a strongly varying complementary torsional moment and its correspondent distortional quantity, which are almost zero anywhere else. The simplified method gives a good approximation of all the derivatives, except a large error on the second derivative of $\theta \left(z\right)$, magnified at the loaded region, consistently with what is noticed in Figure 14.

5.1.1. Stress Analysis in Pure Torsion

A stress analysis was carried out by referring to the Fourier solution with $N=50$ terms. The stresses were evaluated by referring to the uncracked cross section. In order to separate the causes of discrepancies between analytical and FEM analyses, the pure torsion case $\left(\theta \ne 0, \phi \equiv 0\right)$ was considered first, for which the Vlasov beam model is usually assumed to be accurate. Thus, the stresses provided by the analytical model via Equations (24) and (25) were compared with those furnished by a FEM model, in which indeformability of the cross section was enforced via “body constraints”, simulating thin diaphrams at any cross sections of the mesh. It should be remembered, that the analytical model proposed here reproduces the Vlasov model by further encopassing for variation in the stresses within the thickness.

The pattern of stresses on a tipical cross section was examined first. Figure 16 reports, by continuous/dashed lines, the analytical longitudinal normal stresses ${\sigma }_{z,max}^{m+f}:={\sigma }_z^m+{\sigma }_{z,max}^f$ and the tangential stresses ${\tau }_{max}^{m+f}:={\tau }^m+{\tau }_{max}^f$ evaluated at the middle abscissa of the load ($z=7.5 \mathrm{m})$ and at bottom/top edges of the cross section. The bullets represent numerical FEM results. It appears that the “enriched Vlasov model” captures very well the state of stress sufficiently far from the corners of the cross section. Close to them, indeed, and probably exhalted from the abrupt change in thickness, the tangential stresses violate the Bredt law; this occurrence, in turn, entails a local disturbance on the longitudinal stresses for equilibrium reasons. The maximum differences between analytical and numerical results are of the order of 10–15%.

The dependence of the stresses on the longitudinal abscissa was then investigated. Figure 17 shows how the normal stresses at the corners and the tangential stresses at the half side $s_i=0$ depend on z. It appears that the differences between the analytical and numerical models disappear at large distance from the load, thus confirming the local nature of the disturbance.

Concerning the tangential reactive stresses, which are related to the complementary torsional moment, they are negligibly small at $z=7.5 \mathrm{m}$ (not shown here). However, as Figure 18 shows, they are significant at $z=11.25 \mathrm{m}$ (i.e., where the load has a discontinuity point) according to the diagram of ${\theta }^{‴}\left(z\right)$ in Figure 15e.

5.1.2. Stress Analysis in Torsion–Distortion

The complete analytical solution ($\theta \ne 0, \phi \ne 0$) is now considered, and relevant stresses (Equations (26) and (27)) are compared with those furnished by the unconstrained FE model. The membrane longitudinal normal stresses at the corners, ${\sigma }_z^m\left(z,\frac{b}{2}\right)$, ${\sigma }_z^m\left(z,\frac{h}{2}\right)$, are plotted in Figure 19a,b; they are found to be in excellent agreement with FEM analysis. Moreover, they are about three times larger that the normal stresses induced by torsion alone (Figure 17), thus highlighting the role of warping by distortion. The membrane tangential stresses at half side, ${\tau }^m\left(z,s_1=0\right)$, ${\tau }^m\left(z,s_2=0\right)$, are plotted in Figure 19c,d. Here, the contributions of active and reactive parts, together with their sum, are displayed. It is seen that the reactive component associated with distortion is essential, as discussed before, to achieve a good agreement with the numerical results. Finally, the maximum flexural tangential stresses at half side ${\tau }^f\left(z,s_1=0,y_1=\frac{t_1}{2}\right)$ and ${\tau }^f\left(z,s_2=0,y_2=\frac{t_2}{2}\right)$ are plotted in Figure 19e,f, showing, again, an excellent accordance between analytical and numerical results. The flexural component of the tangential stress is found to be comparable or even larger than the membrane component, thus revealing the importance to consider the variation in stresses inside the thickness of the cross section of reinforced concrete bridges.

The transverse normal stresses, mainly related to the distortion, are reported in Figure 20. Stresses acting at the cross section $z=7.5 \mathrm{m}$ are displayed in Figure 20a,b. The linear law predicted by the analytical model is confirmed by FEM to within an error of about 10% at corners. When the variation with z of the stresses at corners was analyzed (Figure 20c,d), a 10% error is confirmed just in the region loaded, while the error almost vanishes far from this zone. As a curiosity, while the analytical model predicts the same absolute values for the stresses at the two faces, the FEM analysis reveals slight differences, more significant on one of the two faces. The reactive transverse normal stresses, which are associated with the normal forces induced by the bending moments, were checked as an attempt to explain this discrepancy, but they were found to be negligibly small. Therefore, such a phenomenon cannot be captured by the unique distortion mode considered here. Moreover, further investigation would be needed to unfold the effects of the way the external loads are applied to the box-girder, which, according to the De Saint Venant principle, are not negligible close to the loaded region.

As a final comment, it is important to notice that the distortional stresses in Figure 19 and Figure 20 are of the same order of magnitude or larger than the stresses due to pure torsion in Figure 16 and Figure 17. Therefore, distortion plays an important role in this sample system.

6. Conclusions

A minimal linear and elastic model of rectangular box-girder, undergoing twist and distortion, was formulated in the framework of GBT. The model is useful for bridge analysis in skew-symmetric loading conditions caused by the eccentricity of live loads. The displacement field was approximated as the superposition of two cross-sectional modes, one describing twist and the other describing distortion, all given in closed form. Modes are modulated by two amplitude functions, defined on the beam axis, governed by a set of two coupled fourth-order ordinary differential equations, which generalize the well-known uncoupled equations by Vlasov for torsion and of beam on elastic soil for distortion. The method, being analytic, is able to supply closed-form expressions for stresses in the whole domain.

Algorithmic aspects related to these equations were tackled, leading to the following conclusions.

• The coefficients of the derivatives of fourth-order are mainly generated by nonuniform warping both in twist and distortion. Therefore, among them, the contribution of the flexural nature of the plates are negligible. In contrast, warping terms, although small, cannot be neglected when the displacement wavelength is short. It is argued that they could be relevant in describing boundary layers, e.g., produced by constraints preventing free warping.
• The distortion-to-twist ratio was proven to be of order 1 for thin and short girders.
• The Fourier analysis is a convenient and efficient tool to analyze simply supported girders warping free at the ends; for other boundary conditions, exact integration of the differential equations should be carried out. By following the literature, a simplified procedure was illustrated, which calls for neglecting all coupling terms, solving two independent problems and superimposing the effects. The two problems are (i) the Vlasov beam under torsion and (ii) the Winkler soil equation-like beam for distortion.

A sample system was considered, for which both Fourier analysis and uncoupled equation integration were carried out. The following conclusions were drawn.

• Fourier analysis works well even for non-smooth loading conditions, provided that a sufficient number of terms is accounted for in the series. The exact integration of the uncoupled equations gives reasonably good results, with errors of about 10% with respect to the coupled Fourier representation.
• Stresses due to torsion mainly consist of (i) normal longitudinal components equilibrating the bimoment; (ii) active tangential stresses, as given by the Bredt theory; and (iii) reactive tangential stresses equilibrating the complementary torsional moment due to warping. All these effects are significant, except for the reactive tangential stresses, of which the influence is appreciable only close to the discontinuity points of the load.
• Stresses due to distortion consist of (i) membrane normal stress in the longitudinal direction, triggered by nonuniform warping, kinematically compatible with the loss of the shape of the cross section; (ii) membrane normal stresses in the transverse direction, generated by the frame-like behavior of the cross-section; (iii) flexural normal stresses in the longitudinal direction, generated by the flexure of the plates associated with the longitudinal modulation of the frame-deflection; and (iv) tangential stresses, generated by the torsion of the plates, made of an active and a reactive component. Among these stresses, (i) and (ii) are the most important. However, among the tangential stresses, the reactive component cannot be neglected.

Finally, all the results have been validated by a Finite Element Analysis. In spite of a remarkably larger computational effort, the differences between numerical and analytical results were found to be small. It is therefore confirmed that the simple analytical model proposed here is able to capture, with acceptable precision, the mechanical behavior of the girder.

The research is susceptible to some extensions. In the static field, (i) the exact solution of the two coupled equation could be found and more complicated boundary conditions could be analyzed and (ii) a one-dimensional finite element could be formulated to account for highly segmented or lumped loads as well for multi-span bridges. In the dynamic field, (iii) the free and forced vibrations of box-girder in torsion–distortion could be analyzed.