Validation of Bredt’s formulas for beams with hollow cross sections by the method of asymptotic splitting for pure torsion and their extension to shear force bending

The equations governing pure torsion of prismatic beams with thin-walled closed cross sections, known as Bredt’s formulas, are verified using the method of asymptotic splitting. In particular, the strong formulation of the Saint-Venant problem of a straight beam is expanded asymptotically. We begin by validating well-known technical assumptions for the shear stress distribution. Furthermore, the influence of a transverse force acting on the beam is considered. This shear force causes a deformation of the cross section, and therefore an adaption of Bredt’s formulas is needed. Two distinct formulations of the shear center, called the kinematic and the energetic shear center, are obtained. The latter is verified in numerical experiments.


Introduction
Due to their lightweight compared to the stiffness, single-and multi-cell thin-walled beams are widely used in engineering. Useful design tools for such structures are the formulas of Bredt, which in the case of a single-cell cross section read In (1), the torque M T is expressed as a function of the shear modulus μ, the torsional rigidity J T , and the rate of twist α. The area enclosed by the center line Γ is defined as A Γ , and the wall thickness h(s) varies along the arc coordinate s. The formulas are derived under the engineering assumption of a constant shear stress over the thickness. As a further assumption, shear stresses orthogonal to Γ are neglected. Equation (2) relates the shear stress τ with the deformation. In [1], the same simplifications are used to obtain the position of the shear center by using Bredt's formulas as an approximation even in the case of combined bending and torsion. We validate the engineering assumptions by an asymptotic analysis of an arbitrary single-cell cross section. Therefore, we consider the Saint-Venant problem of a linear elastic beam. The analytical description of linear elastic beams with arbitrary, also multiply connected cross sections, was done by Eliseev [2,3]. In the limit of small thicknesses h, represented by a formal small parameter λ, the cross section can be characterized by its center line and a small dimension orthogonal to it. The strong form of the equations is developed asymptotically to obtain the principal terms of shear stresses and related quantities. The determination of the principal terms requires the solvability conditions of the minor terms, which is intrinsic to the procedure of asymptotic splitting. This method has also been applied to the actual problem but with constant wall thickness in [3] and to the problem of constrained warping of thin open cross sections as well as plate equations in [4][5][6]. Another possibility of an asymptotic analysis is shown in [7], in which also higher terms are derived, but again only constant thickness is considered. Further asymptotic solutions of the problem were found by the authors of [8,9], who used the weak formulation of the problem.
In the present paper, the solutions are used to verify Bredt's formulas and to find expressions of the shear centers. There are two definitions of the shear center. The first one defines it as the point, in which a single shear force has to act such that no twisting occurs. The second one defines it as the position where the shear force and an additional applied torque do not perform work on each other's deformations. It can be shown that for open thin cross sections these definitions are asymptotically equivalent, whereas in the case of multiply connected ones these two are distinct. In the present paper, we show that even in the case of a single-cell cross section with constant wall thickness the difference remains, a result that so far has not been reported in the literature. The analytical results are validated numerically against the formulation of Lacarbonara [10].

Notations and formulation of the problem
We study stresses and displacements of a Saint-Venant solution for a prismatic beam, whose axis is directed along the unit vector e z with coordinate z, and with the constant cross section placed orthogonal to it in the (x, y)-plane with a position vector x = xe x + ye y (Fig. 1). In case of a thin-walled cross section, we describe this position vector by with the center line x 0 as function of its arc length s, the coordinate n with − h 2 ≤ n ≤ h 2 orthogonal to the center line(with unit vector e n ), and λ being a formal small parameter, representing the smallness of the thickness h.
is the unit tangential vector of the center line Γ of the thin-walled profile. With the curvature κ, we can write the normal vector which is orthogonal to e t . With the Frenet formula the nabla operator is determined by with the Lamé coefficients ∇x = I 2 holds with the identity tensor I 2 . Similar to the vectors e t and e n of the center line, there are a tangential vector t and a normal vector n at the boundary of the cross section. Because of the varying thickness, the directions of t and e t as well as n and e n are different. Note that we have e z = e x ×e y = −e t ×e n . Although the small parameter λ has already been introduced, we present the general exact solution of the Saint-Venant problem. The origin x = 0 is chosen as the area center of the cross section; therefore, with no in-plane components. The equilibrium equations for the axial component σ z and for the shear stresses τ are a being the mean axial stress, and b = ∂b/∂z determines the linear distribution of the axial stress in the cross section. We can write with the tensor of the area moments of inertia J and the shear force Q 0 , which is the in-plane part of Q. Furthermore, we have to satisfy the compatibility condition written in the form of the Beltrami equation, with ν being Poisson's ratio. The shear stress vector can be expressed as with in which μ is the shear modulus and α a constant. Later, we will see that α has the meaning of the average rate of twist of the cross section. The boundary conditions are fulfilled if τ · n = n · ∇ϕ + (e z × n) · ∇ψ = 0.
The first term implies that the normal derivative of ϕ is zero at the boundaries, whereas the second term enforces the tangential derivative of ψ to be zero. From this second expression, we conclude that ψ must be constant at each closed contour of the boundary. In the case of a single-cell cross section, there are two such contours; therefore, we choose one constant to be zero, the other one is yet undetermined.

The method of asymptotic splitting
The asymptotic method used in this paper is based on introducing in the problem a formal small parameter λ, which is first used to obtain a serial expansion for λ → 0, but is set equal to 1 after finishing the procedure (see also in Andrianov [11]). Assuming the order of the principal terms in the solution, we balance the leadingorder terms in the equations. In general, the solution is not unique yet. Therefore, higher orders of λ have to be considered. The method is demonstrated by considering a linear algebraic system which becomes singular, if λ tends to zero. Therefore, we seek u in the form of a power series in λ, starting with λ −1 : Inserting u finds the principal order (λ −1 ) term to The solution u (0) is a linear combination of the fundamental solutions ϕ k , but the coefficients a k are still undetermined. Therefore, we have to consider the first minor order term (λ 0 ): Both sides are multiplied by ψ T i , and the solutions of the conjugate system C T 0 ψ i = 0; hence, we obtain This further condition results in a linear system for the coefficients a k . If this system is again singular, further minor terms have to be considered. For a more detailed presentation of regular and singular perturbed linear systems, see [12]. Depending on the problem, also boundary layers may arise [13]. Note that instead of introducing λ it would be possible to nondimensionalize the equations and search for a dimensionless parameter, which is physically small. But in our opinion the introduction of λ is, at least when searching just the principal terms, more convenient.

Bredt's formulas of a single-cell cross section
First, the shear stress is determined to validate that it is constant over the thickness, and next the displacements are used to identify this constant.

Validation of the assumptions of the shear stress τ
The equilibrium equations and the compatibility conditions are equivalent to two partial differential equations (15), (16), namely Poisson equations for the functions ϕ and ψ. From (12), we conclude with Seeking the solution for ϕ as and inserting into (15), we find With a series expansion in terms of λ, this equation can be arranged by orders of λ −3 , λ −2 and λ −1 . Hence, we obtain The boundary conditions on the outer and inner surface are n · ∇ϕ n=±h/2 = 0. With the boundary condition becomes Collecting terms of different orders in the boundary condition yields λ 0 : ± ∂ϕ (2) ∂n Together with (27), (28), and (29), we conclude that ϕ (0) = ϕ (0) (s) and ϕ (1) = ϕ (1) (s), as well as Inserting ϕ (2) into the two boundary conditions of order λ 0 at n = ±h/2 gives − 1 2 ∂h ∂s such that C 1 (s) must vanish, and we have from which ϕ (0) can be computed. To start the asymptotics of ψ, it is imperative to know the order of the right-hand side terms. Therefore, we compute the torque to Knowing the result for the special case of pure torsion of a circular ring , we conclude from (38) that the principal order of ψ has to be λ 0 and the rate of twist has to be α = λ −1 α (0) + · · · . Therefore, the boundary conditions are ψ(n = +h/2) = 0 and ψ(n = −h/2) = λ 0 C (0) + · · · . With ψ = λ 0 ψ (0) + · · · , the asymptotic equations are with the principal solution Finally, one can insert the results into (14) to compute the shear stress vector Hence, the principal term of the shear stress vector, is constant through the thickness and tangential to the center line. With b = 0 in pure torsion, ϕ (0) = 0 holds. Setting λ = 1, the shear stress vector is τ h = C (0) e t , and C (0) is computed from the circulatory theorem With the shear flow T = τ h = C (0) Bredt's formulas are obtained. In this paper, we are further interested in the extension of Bredt's formulas to account for shear forces. With the definition of the shear flow the integration of τ across the thickness yields Inserting this formula into (37) yields the known equations ∂ T ∂s with b (0) · x 0 = ∂σ z /∂z , from which we find the shear flow Next, we establish a further condition to evaluate T 0 = T (s = 0), which follows from the continuity of the displacements.

Consideration of the displacements
Due to the absence of in-plane stresses in the cross-sectional plane, the in-plane strains ε 0 = (∇u 0 ) sym = −νI 2 ε z are caused only by the Poisson effect, with the 2D identity tensor I 2 and the axial strain ε z . The resulting displacements u = u 0 + u z e z according to Hooke's law applied to the principal strains and the axial stress σ z , including an additional in-plane rigid body motion U 0 + ωe z × x as well as an additional axial displacement U z (x), are For details of this classical result, see [2]. The displacement U z (x) may be used, too, to obtain the solution, but as seen below, usage of the stresses U z (x) is not needed. The displacements are inserted into Hooke's law, If we insert the displacements as well as ψ and ϕ into Hooke's law and take the rotor ∇×, we find α = ω ; see Eliseev [2]. This is the proof of the above statement that α is the mean rate of twist. It should be mentioned that because of Poisson's ratio the rate of twist is not constant within the cross section. The resulting equation splits into a part, which depends on z and one part, which is independent of z. Both must hold on their own. The part dependent of z is In the following, we consider only the part of the equation which is independent of z, which reads Next we take the principal terms of the equation and integrate along the center line Γ , with A Γ being the area enclosed by the center line. With ∂x 0 /∂s = e t and integration by parts, we find where we already can see (2) in the case of pure torsion (Q 0 = 0). According to the Stokes theorem, the integral on the right-hand side can be transformed to an area moment of first order, which does not vanish, because the coordinate system was set in the gravity center of the cross section, which is not the gravity center of the embedded area. This leads finally to an expression for the remaining unknown C (0) , with A Γ x Γ = A Γ xd A Γ . According to Eliseev [3], the work conjugate of the moment M z is Inserting θ z into the above equation results in Further, by inserting τ = T / h, we obtain T 0 as With (48), we can write the leading order λ 0 of the torque as with p = e n · x 0 . Integrating by parts and noting x 0 hds = 0, we obtain With the abbreviations The terms in the brackets represent the principal term of the warping function W from pure torsion; therefore, we can write If we assume pure torsion (b = 0), the rate of twist α (0) is equal to θ z , and we get which validates (1). Repeating (61), the adapted second formula of Bredt including shear forces is In the case of pure torsion, with θ z = α (0) , we recover (2).

Shear center
The kinematic shear center x * is defined as the point where a shear force Q 0 causes no twisting of the beam (α = 0). The torque resulting from the shear stresses τ has to compensate the torque caused by the shear force, that is, Next, we substitute M z considering α = 0 and get We substitute b (0) from (12), cancel Q 0 , and multiply the resulting equation with ×e z . The resulting equation is Expressions for the coordinates of the shear center can be written in a simple form if the principal axes of inertia coincide with the axes of the coordinate system, because J = J y e x e x + J x e y e y , and therefore also its inverse is diagonal. For the x-component, we then write In contrast to the kinematic shear center, the energetic shear center x * * is defined at the point where an acting shear force has no influence in the work conjugate of the torque M z , which implies that θ z has to be zero. This results in respectively, the x-coordinate The two shear centers only differ in the term containing ν. They coincide if x Γ = x A (= 0). The difference also disappears if ν = 0, which represents rigidity of the projection of the cross section onto the cross-sectional plane.

Numerical example
Next, the above theoretical considerations are validated numerically. We consider two cross sections symmetric to the x-axis (Fig. 2). The first one is a thin hollow square profile with the right wall being βh thick and the other three walls being h thick. The second one has the shape of a hollow C-profile with a constant wall thickness. For the first one, we study the convergence of the asymptotic solution; the second example was chosen to proof the statement that for thin single-cell cross sections with constant thickness the difference of the two shear centers vanishes is wrong [14]. The computation is done by the above analytical formulas as well as numerically. We approach the exact three-dimensional solution by using the method of finite differences. The numerical solution makes no assumption about the thickness of the cross section. Instead of solving the above Laplace equations for the stress functions (15) and (16), the basis of the finite difference scheme are the equations established by Lacarbonara [10] (only the x-component is needed in the two symmetric examples): This approach is called the displacement approach, because it is based on computing the divergence of both sides of (55) and representing the displacement U z (x) as a superposition of three warping functions W 1 , W 2 , and W : It can be seen easily that W is the well-known "ordinary" warping function from the problem of pure torsion. It should be noted that the above equations for W 2 are only valid if the coordinate axes are principal axes of inertia. The advantage of the displacement approach of Lacarbonara compared with equations (15) and (16) is that the boundary condition is a Neumann one, which is known at all boundaries, whereas the above formulation uses the Prandtl stress function, which is constant at the different boundaries, but a priori unknown and has to be calculated with compatibility conditions (see above for C (0) ).

Finite difference scheme
The numerical solution is obtained by discretizing the Laplace equation in the two-dimensional domains (Fig.  2) using the finite difference method. To keep it simple, the differences Δx and Δy are constant. In the case of a node (i, j) inside the domain, the scheme requires five points and reads At the boundaries, the scheme has to be adapted. In our case, the normal vectors are only horizontal or vertical which simplifies this process. Therefore, we have four different kinds of boundary conditions, depending on the direction of the normal vector. For example, the boundary conditions at a straight boundary with constant x = x i BC and n y = 0 are Analogous relations are obtained for the other boundaries as well. Hence, the numerical scheme is complete. At a certain boundary, the appropriate node is substituted by one of the boundary conditions. If the middle node (i, j) is a convex corner, two points have to be replaced. Concave corners are treated like inner nodes.
Because there are only Neumann boundary conditions, the solution is determined up to an additive constant, and the solution for an arbitrary point is set to zero, e.g., the left-low corner W 0,0 2 = 0. Note that in (75) the term vanishes, which confirms the arbitrariness of W 0,0 2 .

Numerical results
The solution is obtained for parameters ν = 0.3, a = 1, β = 2 and different thickness values h. We choose Δx = Δy and at least six points per wall thickness of the thinner walls. Then, the convergence is tested by repeated reduction of Δx and Δy by the factor 2. As we can see in Fig. 3, for small wall thicknesses (large ratios a/ h), the two centers, measured from the left outer edge, remain distinct, and the numerical solutions tend toward the asymptotic expressions. Note that x * and x * * are computed from (70) and (72), added by x A , measured from the left edge. Next, the error due to variable thickness of the cross section is discussed. Therefore, the energetic shear center x * * of the quadratic cross section is presented as a function of the wall thickness ratio β. Three thickness values h = a 500 , a 100 , a 20 , with a = 1 are considered. In Fig.4, the asymptotic results are compared with the numeric ones. For β = 1, the cross section is double-symmetric, and therefore both solutions are independent of h exactly at the center of the cross section, a fact which would count for the kinematic shear center, too. With increasing β, the shear center at first shifts to the right up to a maximum value and decreases again. The maximum value of the shear center of the curve h/a = 0.05 is at a point at which the thickness of the wall is about 25% of the square's length.
After this numerical validation of the asymptotic equations for the shear centers, we proceed with a further instructive example, showing that also for a warping-free cross section x Γ = x A is possible. It is well known that cross sections with constant thickness, which represent tangential polygons, are warping-free. Therefore, .
We can see that only in the case of γ = π 3 , which represents an equilateral triangle, the two expressions are equal.

Conclusions
In this paper, we studied the asymptotic behavior of the solutions of the Saint-Venant problem for a linear elastic beam with a hollow cross section under the action of torsion and shear force bending. Assuming thin-walled cross sections, we introduced the small parameter λ, representing the thinness of the walls. After expressing the strong form of the equations in terms of λ and merging terms of equal order, we first managed to validate the engineering assumptions. Further, we got well-known formulas, like for the shear flow and finally Bredt's formulas for the torsion of rods. We noticed that if the Poisson effect is not neglected within the cross section, the second formula of Bredt has to be adapted by using the energetic twist rate θ z instead of the kinematic one. The leading-order terms of the kinematic and energetic shear centers were derived assuming a variable thickness of the cross section. In the final numerical investigation of two hollow cross sections, the exact location of the shear centers was calculated for different thickness values. A good agreement with the obtained asymptotic formulas was observed when the thickness is small. Furthermore, it has been shown that the energetic and kinematic shear centers are in general distinct even in the case of a single-cell hollow cross section with a constant small thickness. The opposite statement, which is sometimes met in the literature [14], is thus proven to be incorrect.