On incorporating warping effects due to transverse shear and torsion into the theories of straight elastic bars

By endowing El Fatmi’s theories of bars with first-order warping functions due to torsion and shear, a family of theories of bars, of various applicability ranges, is effectively constructed. The theories thus formed concern bars of arbitrary cross-sections; they are reformulations of the mentioned theories by El Fatmi and theories by Kim and Kim, Librescu and Song, Vlasov and Timoshenko. The Vlasov-like theory thus developed is capable of describing the torsional buckling and lateral buckling phenomena of bars of both solid and thin-walled cross-sections, which reflects the non-trivial correspondence, noted by Wagner and Gruttmann, between the torsional St.Venant’s warping function and the contour-wise defined warping functions proposed by Vlasov. Moreover, the present paper delivers an explicit construction of the constitutive equations of Timoshenko’s theory; the equations linking transverse forces with the measures of transverse shear turn out to be coupled for all bars of asymmetric cross-sections. The modeling is hierarchical: the warping functions are numerically constructed by solving the three underlying 2D scalar elliptic problems, providing the effective characteristics for the 1D models of bars. The 2D and 1D problems are indissolubly bonded, thus forming a unified scientific tool, deeply rooted in the hitherto existing knowledge on elasticity of elastic straight bars.


Introduction
The paper concerns modeling the spatial deformation of straight and prismatic elastic bars, of both solid and thin-walled cross-sections, including warping phenomena due to transverse shear and torsion. Introducing the warping effects necessitates refining the kinematical assumptions to encompass all the deformation modes due to: tension, bending in two orthogonal planes, transverse shear in two directions, and torsion. The process of refining the kinematical assumptions concerns mainly the axial displacement field, the transverse displacements being assumed like in Vlasov's [1] theory; they describe the rigid motion of the cross-section: its transverse displacement and rotation around the bar axis. A variationally consistent theory of bars including all the deformation modes mentioned above and based on the properly constructed kinematical assumptions has been proposed by El Fatmi [2,3]. In the present paper, we refer to this approach in both its versions: general (FG) and simplified (FS). According to the latter approach, the warping is determined by the measures of transverse shear and torsion. Within this theory, the primal (kinematic) unknowns are the same as in Timoshenko's theory: three displacements and three angles of rotation. The theoretical schemes of both the FG and FS theories involve neither a construction of the center of shear (coinciding with the center of torsion if Poisson's ratio is set to zero) nor a construction of the warping functions; the general scheme of the theory may involve any forms of the warping functions. The present paper is aimed at specifying this theory by assuming the warping functions of the first-order approximation due to shear and the commonly accepted pure torsion warping function. Under this choice of the warping functions, the constitutive equations linking the transverse forces and Timoshenko's measures of the transverse shear turn out to be coupled, as noted in the paper by Mason and Herrmann [4], corrected by Schramm et al. [5], and then discussed in Dong et al. [6][7][8] and Romano et al. [9], while the warping functions implicitly assumed by El Fatmi [2] result in a decoupling of these equations, see Eq. (37) therein.
In contrast to the warping due to shear, the description of the warping due to torsion is commonly accepted as that known from the theory of pure torsion, see Sect. 2.2. in Gruttmann et al. [10]. The functions of warping due to shear may be constructed in several manners, see Gruttmann et al. [10] and Gruttmann and Wagner [11]. In the present paper, their simplest version will be assumed, based only on the equilibrium equations, disregarding the compatibility equations. They will be called first-order approximation warping functions; they correspond to the second-order approximation upon putting zero value of Poisson's ratio. To construct the warping functions, the numerical codes in C++ have been prepared and implemented. To make the present paper self-contained, the numerical approach is thoroughly presented (Sect. 4.4 and Appendix) to clear up in detail how to solve the specific numerical problems occurring while solving the auxiliary scalar elliptic problems, the solutions of which precede the construction of the warping functions. Let us stress that the whole approach is hierarchical, it consists of three steps: (a) constructing the necessary auxiliary fields within the 2D domain (i.e., defined in the cross-section of the bar) with appropriate accuracy and compute the necessary cross-sectional characteristics, (b) solving the static problems within the 1D theory of bars, (c) recovery of stress and strains within the bar viewed as a 3D body. The present paper is focused on the stage (a) and provides the 1D models. The methods of solving the latter and the stage (c) are not subject of the present paper.
The El Fatmi model may be improved by incorporating more accurate warping functions, as proposed in Dikaros et al. [12], Sapountzakis and Mokos [13], Sapountzakis and Dourakopoulos [14], yet such more involved models will not be discussed. One of the reason is that by virtue of neglecting Poisson's effects in the modeling of the warping due to shear the two points: shear center and torsion center become one point; indeed, these two points coincide only if the mentioned Poisson effect is disregarded, see Romano et al. [9]. Considering this effect would mean a great price to pay; the consequences would be far-reaching; practically, the whole simplicity of the modeling would be lost. Let us stress: the main aim of the present paper is to construct the theories of bars and not merely refine computing the stiffnesses, which can be done a posteriori. The warping functions used here enter the kinematic assumptions at the very beginning, hence their maximum simplicity is highly advisable.
By neglecting the warping due to shear in the FS model one comes across the Kim and Kim [15] theory. However, in its original version this approach is not consequently displacement based; it is constructed by using the stress-displacement (saddle point) functional; the construction of stiffnesses is based on the Vlasov formula for stresses within the thin-walled bars. Since the aim of the present paper is to construct a theory applicable for bars of arbitrary (both solid and thin-walled) cross-sections, the original Kim and Kim method cannot be used here. One of the exceptional properties of the Kim and Kim-like theory proposed is its displacement-based modeling without any reference to Vlasov's approximations pertaining to the thin-walled profiles.
El Fatmi's [2] general theory is a convenient basis for constructing the theory of Librescu and Song [16], see Sect. 3 therein, which is an extension of Vlasov's theory towards including transverse shear deformations, see also Pavazza [17]. This theory is discussed in Sect. 8.
Alternative theories of bars can be found in the papers cited in El Fatmi [2] and in the recent papers discussing specific problems of non-prismatic and composite bars, see Schneider and Kienzler [18], Endo [19].
As mentioned above, the present paper delivers a variational derivation of the constitutive equations linking the shear forces T z , T y with the measures of transverse shear γ z , γ y , the coupling (cross-property) stiffness being derived consistently (Sect. 5.5), also disclosing (and not imposing) symmetry of these equations. The derived stiffnesses due to transverse shear occur to coincide with those found by Schramm et al. [5], see equations (51-54) therein, upon putting ν = 0. However, the paper by Schramm et al. [5] does not deal with torsion, while the present paper delivers the proof that the constitutive equations for the transverse forces and the torsional moment are decoupled within the FS approach augmented with the first-order warping functions.
The coupling of the constitutive equations for the transverse shear forces shows that for each test of bending a lateral displacement must occur accompanying the bending in the plane orthogonal to the plane of the load. In particular, we conclude that the theories of thin bars (like Bernoulli-Euler's and Vlasov's) which neglect the transverse shear deformation are incapable of describing this phenomenon. This suggests that the mentioned thin bars theories should not be used if the cross-sections of bars are asymmetric. However, the analysis in Sect. 13.1 shows that this effect is negligible, at least within the linearized setting.
The paper is concerned with the straight and prismatic bars, hence of cylindrical shapes, of the axis x and of the cross-section domain A (its area will be denoted by the same font, which should not lead to misunderstandings). Thus, the vectors n outward normal to the surface of the bar are orthogonal to the axis x. The domain A may be multi-connected. It will be parameterized by the Cartesian orthogonal system (y, z); the system (x, y, z) is assumed as right-handed. The components of displacements (u x , u y , u z ) determine the components of strain, where (·) ,x = ∂/∂ x, (·) ,y = ∂/∂y, (·) ,z = ∂/∂z. According to the kinematic assumptions imposed later, the remaining strain components: ε y , γ yz , ε z , vanish. The material of the bar is assumed to be isotropic and homogeneous. The constitutive equations will be assumed in the form commonly accepted in the papers on the elasticity of thin and straight bars: where E represents Young's modulus, E = 2(1 + ν)G, G is the shear modulus, while ν is Poisson's ratio. The body forces are omitted. The first equilibrium equation reads The following conventions will be adopted. The scalar product of two fields f and g given on A is defined by where d A = dydz. The matrix (y|y) (y|z) (z|y) (z|z) of moments of inertia is always positive definite, which follows from the Cauchy-Schwarz inequality. The following notation is used for brevity: The axes y, z will be chosen as principal axes of the domain A, or The point O of coordinates (0,0) is the center of gravity of A, while the principal moments of inertia are: J y = (z|z) , J z = (y|y). The cross-sections may be: bisymmetric (both axes y and z are axes of symmetry), mono-symmetric (the cross-section has a single axis of symmetry which is also one of the principal axes) or asymmetric (neither of the principal axes y, z is a symmetry axis). The contour of the plane domain A will be denoted by ; the unit vector n within this plane, outward normal to the contour has components n y , n z . The directional derivative of the function f along the vector n is defined by ∂ f ∂n = f ,y n y + f ,z n z .

The center of torsion and the warping due to pure torsion
To make the paper self-contained, let us recall the definition of the warping function accompanying the pure torsion. We follow the paper by Gruttmann et al. [10]. In the first step, we construct the function ω o = ω o (y, z) being the solution to the following elliptic problem: The solution exists and is unique under well-known regularity conditions concerning the domain A. Let us define a new function: where (y S , z S ) are coordinates of a certain point S. Let us introduce the requirements Due to the assumption (y|z) = 0 of y, z being the principal inertia axes, we find The point S will be called the center of torsion. In case of neglecting Poisson's ratio effect on the form of the warping due to shear, the center of torsion coincides with the center of shear, see Gruttmann et al. [10] and Romano et al. [9]. The axis x links the centers of gravity of the cross-sections A forming the axis of bending. The centers of torsion/shear S form the axis x S , being the axis of torsion (or, equivalently, the axis of transverse shear). The function (2.2) satisfies the variational equation The auxiliary functions modeling the warping due to transverse shear The warping due to shear will be modeled by the warping functions of first order, neglecting Poisson's ratio effect, see Gruttmann et al. [10] and Gruttmann and Wagner [11]. Let the functionζ =ζ (y, z) be the solution to the elliptic problem: If the last condition in (3.1) is neglected, its variational form reads where ∇ζ · ∇v =ζ ,y v ,y +ζ ,z v ,z . Let the functionη =η(y, z) be the solution to the elliptic problem: If the last condition in (3.3) is neglected, its variational form reads.
Substitution of v =η into (3.2) and v =ζ into (3.4) results in the Betti-like identity of reciprocity, Let us define we conclude, by using the Cauchy-Schwarz inequality, that D>0. Let us introduce the new fields defined on the domain A. As will be made clear later, the functions η(y, z) − y, ζ(y, z) − z model the warping due to transverse shear in the directions y and z.

Construction of the warping functions
Let us introduce the first-order warping functions The function g 1 (y, z) models the warping accompanying the torsion induced by a constant torque. The functions g 2 (y, z), g 3 (y, z) are the warping functions due to shear in the z and y directions, induced by the given transverse forces at the ends of the bar. The functions (4.1) are orthogonal to the functions: 1, y, z with respect to the scalar product (1.4) These functions cannot be formed simpler, having in mind the equilibrium equations to be satisfied for the static problems of homogeneous states of torsion and shear induced by the constant loads at the ends of the bar. The functions (4.1) neglect the Poisson's ratio effect. Due to this simplification, the two points: center of shear and center of torsion coincide (point S), which essentially simplifies the construction of the theories of bars. Two (and not three) points: O and S in the domain A will be singled out. Thus, the formation of a correct theory of bars must be preceded by solving the three elliptic problems (2.1), (3.1), (3.3) on the plane domain A. The numerical construction of all the warping functions will be presented in Sect. 4.4.

Warping accompanying the pure axially uniform torsion
Consider the state of displacements where ρ = dθ(x)/dx. Assume now that the torsion is axially uniform, or According to (1.2), the associated stresses are: Hence τ xy,y = Gω ,yy ρ o , τ xz,z = Gω ,zz ρ o . Thus, the equilibrium equation (1.3) is satisfied if ω = 0, which holds, see (2.1), (2.2). If the surface of the bar is free of any load, then σ x n x + τ xy n y + τ xz n z = 0. (4.5) The bar is prismatic, n x = 0, hence the function ω satisfies the boundary condition of the form ω ,y n y + ω ,z n z + (−z + z S )n y + (y − y S )n z = 0, where represent the measures of transverse deformations, known from the Timoshenko bar theory. Let us assume that the transverse shear forces are constant along the bar axis; then the shear deformations are also constant, and then: γ z = γ o z , γ y = γ o y , the quantities γ o z , γ o y being constant. Such a displacement field is associated with stress fields of the form (4.10) and the remaining stress components are assumed as negligible. Like in the Timoshenko theory, we assume that the shear forces are linked with the bending moments by Thus, the bending moments are viewed as being linearly varying along the bar's axis. The constitutive equations are assumed in the form (4.12) where A z , A zy , A yz , A y will be fixed such that the equilibrium equation (1.3) is fulfilled identically. By inserting κ y (x) = M y /E J y , κ z (x) = M z /E J z into (4.9) and then into (1.3), we get (4.14) The equations above deliver four equations for three unknowns, since the relation A zy = A yz should hold. It turns out, however, that we proceed correctly and no contradiction will appear. We substitute (3.9) into (4.14), and, by using (3.1), (3.3) we arrive at In virtue of the identity (3.5), we obtain A zy = A yz , which proves the correctness of the derivation. Moreover, we note that the matrix of the effective areas of the cross-section is positive definite, since: The formulae (4.15) for the effective areas A z , A y , A yz , A zy coincide with the formulae (51-54) by Schramm et al. [5] in case of Poisson's ratio set to zero. The shear correction factors are defined by: k y = A y /A, k z = A z /A, k yz = A yz /A. Impressed by Romano et al. [9], we put forward the conjecture that the spectral radius of the matrix k = A −1 A is not larger than 1. The matrix inverse to A has the form: , (4.17) its form being simpler than (4.15). We conclude that if γ z , γ y are constant along the bar's axis and if the constitutive equations have the form (4.13), (4.15), then to the state of displacement (4.7) a state of stress is associated such that the components of stress: σ x , τ xy , τ xz satisfy the equilibrium equation (1.3) identically. Let us note that the frequently used assumption A zy = 0, A yz = 0 would be wrong, it would lead to mutually contradictory conditions, and not all of the conditions (4.2) and (1.3) could have been satisfied identically. This problem has been, e.g., emphasized in Dong et al. [6].

The state of axial deformation
The state of displacement

Numerical construction of the warping functions
The aim of this Section is to give details of the numerical construction of the solutions to the elliptic problems (2.1), (3.1), (3.3) and of the warping functions formulated in their terms. For the purpose of clearing up the numerical method, we introduce here a specific notation. First, we write: y = (y 0 , y 1 ) = (y, z) ∈ A, ∂u ∂n = ∂u ∂ y 0 n 0 + ∂u ∂ y 1 n 1 = ∇u · n, n = (n 0 , n 1 ) = n y , n z ∈ R 2 , (4. 19) and now the three scalar elliptic problems mentioned above can be rewritten by a single scheme, for problems (2.1), (3.1), (3.3), respectively. Let V be the space of test functions defined on A whose first derivatives are piece-wise continuous (the rigorous mathematical definition of V will not be recalled). Let us introduce the bilinear form a (·, ·) : V ×V → R and the linear form L (·) : V → R by (4.22) The variational formulation of (4.20) reads: find u ∈ V such that u = 0 and and if A is a polygon, it may be replaced by the approximate problem: Here V h ⊂ V , being a subspace of V , is spanned finite element-wise by the polynomials of first order, ∀y ∈ R 2 p = p (y) = p 00 + p 10 y 0 + p 01 y 1 , p 00 , p 10 , p 01 ∈ R. where u e i are the unknown values of the function u e h (·) at three subsequent vertices of the triangle A e , while the first-order polynomials N e i : A e → R, i = 0, 1, 2 are the shape functions which depend explicitly on the Cartesian co-ordinates z e i = z e i0 , z e i1 ∈ R 2 , i = 0, 1, 2, of the three vertices defining a triangular finite element A e , see Fig. 1. In the case considered, the formulae defining the shape functions are relatively simple. However, even here it is thought appropriate to avoid using the functions N e i (·) in (4.26) and replace them by (by far more simple) shape functions: defined on the master element; in our problem these functions are expressed by Implementation of the shape functions (4.28) for an arbitrary element necessitates introducing a bijective family of mappings, which link the master element A with an arbitrary element A e such that where (see Fig. 1) The formulae (4.30), (4.31) make it possible to replace Eq. (4.26) by a much simpler one, The definition of the mapping F e = F e 0 , F e 1 is given in the "Appendix", where also other hints of how to implement the finite element algorithm for solving the problem (4.24) can be found.
Let n be the number of the degrees of freedom, here equal to the number of nodes of the FEM mesh. By applying the standard FE procedure to the variational problem (4.24), one finds the linear system K u h = L, (4.33) in which vector L has the dimension n, while the square matrix K of dimensions n × n is singular, since the additional condition u * h = 0 has not been used. To make the solution u * h ∈ R n unique, we apply the above condition thus modifying L to the n+1-dimensional vector L (in which the last component is zero) and modifying K to the form of the rectangular matrixK of dimensions (n + 1) × n. The unique solution u * h ∈ R n of the overdetermined linear system is constructed upon implementing the algorithm SVD (singular value decomposition) in the second author's FEM program. Having u * h along with the interpolation formula (4.32) makes it possible to compute the values of the scalar field u * h at arbitrary point y of the finite element A e . In particular, one can perform the numerical integration with using this field to compute the integral characteristics of the bar's cross-section.

El Fatmi's simplified theory equipped with the warping functions of first-order approximation
Vlasov's [1] theory of thin-walled bars involves the torsional warping effects and neglects both: the transverse shear and the warping due to shear. El Fatmi [2,3] extended this theory to encompass these phenomena. In the general version of the El Fatmi theory, the mentioned three warping deformations are modeled by three independent kinematic variables; in particular, the variation of the warping due to torsion is measured by the field ϑ(x), while in Vlasov's theory this field is equal to the measure of torsion ρ( being the angle of the axial rotation. In the simplified El Fatmi theory, all the measures of warping are exactly equal to the deformations: due to torsion and transverse shears in the y, z directions. The aim of the present section is to construct the simplified El Fatmi theory by equipping it with the first-order warping functions (4.1).

The kinematic assumptions and the associated strains and stresses
In general, the bars undergo the transverse shear, bending, tension/compression and torsion. If the bar is loaded by the forces: N , T y , T z , and by the torque M at its ends, then the displacement fields (4.3), (4.7), (4.18) are accompanied by the stresses satisfying the equilibrium equation (1.3) pointwise. Now we shall assume that the state of displacement has the same form in a general case; we assume the displacement state being the superposition of the states (4.3), (4.7), (4.17), or where the functions ω, ζ, η are constructed by solving the elliptic problems (2.1), (3.1), (3.3), see (2.2), (3.9). Let us note that the orthogonality conditions (4.2) deliver the geometric interpretations of the quantities u, ϕ, β, v, w. Indeed, orthogonalization of (5.1) by the functions: 1, y, z gives u( is an average axial displacement; it is not equal to u x (x, 0, 0), but it will be convenient to treat u as such, since the conjugated end-force N will be treated as being directed along the axis x linking the gravity centers of the cross-sections. Moreover, the quantities ϕ, β are averaged angles of rotations around the axes: −z and y; such definitions are well known from Reissner's theories of plates of moderate thickness, see Eq. (12) in Reissner's paper [20] and Sect. 5.7.1. in [21]. The transverse displacements are simply interpreted as displacements of the shear/torsion center S, since u y (x, y S , Let us recall now the formula for the angle of small rigid rotation around the axis orthogonal to the plane (y, z) : zy = (u z,y − u y,z )/2. Substitution of (5.1) gives zy = θ(x), which confirms that θ(x) represents the angle of rotation of an arbitrary point around the bar's axis. From the theory of plane rigid motion, it is known that the angle of rotation is independent of the choice of the pole. It is convenient, however, to consider θ(x) as referred to the axis x (s) linking shear centers (called the axis of torsion), see Fig. 2. Within the theory considered the following fields play the role of strain measures: Like in Bernoulli-Euler's theory of bars, the quantity ε represents the axial deformation; γ z , γ y are measures of transverse shear; κ y , κ z stand for the measures of flexure; ρ is the measure of torsion; κ ω describes the warping due to torsion; γ z , γ y represent warping due to transverse shear. The notation κ ω is taken from Vlasov's [1] theory of bars.
The strain components associated with the displacement field (5.1) are expressed by The stresses are given by (1.2).

The internal and end equilibrium equations of the bar
Let an overbar indicate a virtual (or test) quantity in the variational equilibrium equation: The loading acting along the surface of the bar determines the virtual work of the span load. This virtual work will be represented by a simplified formula as follows: where p represents the intensity of the axial load reduced to the axis x; q z , q y stand for the intensities of the transverse loads reduced to the axis x (s) ; m is the intensity of the torsional load reduced to the same axis. Assume that the section x = 0 is fully clamped and the end x = l is free, loaded by the tractions of intensities t x (y, z), t y (y, z), t z (y, z). The virtual work of the tractions reads: Substitution of the kinematic assumptions (5.1) rearranges this virtual work to the form involving the end-forces and end-moments expressed bŷ Let us stress thatM y ,M z represent the external moments directed along the axes z and −y, whileM acts along the axis of torsion x (s) ; the lines of actions of the transverse forcesT z ,T y go through the axis x (s) ; their application point is the torsion/shear center S.
The terms composing the virtual work (5.7) can be grouped in an other way, but the expression (5.7) seems to be convenient. Let us define the stress resultants Note that the quantities N , M y , M z refer to the axis x, while the quantities T z , T y , M refer to the torsion axis x (s) , see Fig. 2.
The virtual work of stresses (1.2) entering the l.h.s. of (5.4) equals Substitution of the expressions for the strains (these expressions refer to the virtual strains and unknown strains as well) results in Performing integration over the domain A with using the definitions (5.9) leads tō Now we are ready to construct the equilibrium equations of the theory, along with the natural boundary conditions at the end x = l. By equating the virtual works:L int andL span +L ends , and taking into account arbitrariness of the virtual kinematic fieldsū,β,w,φ,v,θ within the interval (0, l), their vanishing at x = 0 and their arbitrariness at x = l we obtain:

The natural boundary conditions at
where The variational equilibrium equations (5.13) (where the virtual work of end forces has been omitted for simplicity) determine the local equations of equilibrium. The local form of these equations will not be given here; they are easy to derive. The same equations determine the admissible natural boundary conditions of the theory: the conjugated pairs of end-forces and the corresponding displacements are: Thus, at the end of the bar given are: either N or u,…, either B y or ϕ + dv dx . The number of admissible types of the boundary conditions is 2 9 . We understand that also other types of the boundary conditions are possible, if, e.g., the supports are rotated with respect to the axes: y, z. The corresponding admissible end-forces can be found by appropriate rearrangement of the virtual work (5.7).

The constitutive equations
Substitution of the stresses 1.2 along with the strains corresponding to (5.1) into (5.9) gives the following formulae linking the stress resultants (5.9) with the deformation measures (5.2): where S i j = (g i |g j ), and S 11 = (ω|ω) is usually denoted by I ω and called the (torsional) warping moment of inertia. Moreover, the stiffness due to torsion GJ is determined by The identity (2.6) reduces the above formula to the form We see that the torsional stiffness GJ is expressed by the values of the warping function ω along the contour .
We shall prove now that: By (2.5), the first term of the right-hand side of the expression above (or the coefficient at ρ) vanishes. Substitution of (3.9) into the remaining terms of this formula gives, subsequently, the terms G A yz γ z , G A y γ y , where A yz , A y are given by (4.15). In an analogous manner way, we derive the equation for T z . Substitution of the stress components determined by the strains (5.3.1) into the definition of the torsional moment (Eq. (5.9.6)) and using the variational equations of Sect. 3 reduces the constitutive equation for the torsional moment to the form depending only on the torsional measure ρ or M = G Jρ.
The matrix S appearing in (5.17) has the mathematical structure of a Gram matrix and hence is positive definite. We know already that the matrix A is positive definite, see Sect. 4. Thus, the density of elastic energy given by is a positive definite quadratic form of the strain components (5.2).

Analysis of the equilibrium equations
Substitution of the constitutive equations (5.17) into (5.12) determines the bilinear form of the static problem considered. Let us define the collections of the displacement fields: unknown u = u β w ϕ v θ and virtual: u = ūβwφvθ . The virtual work of the internal forces can be recast as a bilinear form: The right-hand side of (5.4) determines a linear form f (ū). We write down Eq. (5.4) in the standard form for eachū kinematically admissible (5.24) referring to the celebrated Riesz' representation theorem, as the bilinear form (5.23) represents a scalar product in appropriately defined space of admissible displacements. Indeed, a F S (u, v) is symmetric and nonnegative: a F S (ū,ū) ≥ 0, since W ≥ 0, and the equality holds if . . , c 6 being arbitrary constants. These rigid body motions describe translations and small rotations; they will be eliminated if the bar is correctly supported. Then a F S (ū,ū) > 0 ifū satisfies correct kinematic boundary conditions. These conditions suffice to prove uniqueness of the solution of the problem (5.24). Moreover, the forms a F S (u,ū) , f (ū) are continuous for the natural regularity conditions imposed on the loading. Consequently, by Riesz' representation theorem, the solution of (5.24) exists; the details of the proof exceed the framework of the present paper.  (−y, z). Thus, c = z|η = y|ζ = 0 and D = y|η z|ζ , hence It is seen that g 2 (y, z) = g 2 (−y, z) , g 3 (y, z) = −g 3 (−y, z), and g 1 (y, z) = −g 1 (−y, z) since now g 1 = ω o (y, z) − z S y. According to (4.15) A yz = A zy = 0, A y = (J z ) 2 /(y|η), A y = (J y ) 2 /(z|ζ ). Consider now the matrix S, S i j = (g i |g j ); now S 13 = 0 since both the functions g 1 , g 3 are odd with respect to y. Thus, according to (5.10) ⎡ We see that a coupling between the torsion and shear along the axis y remains. The other constitutive equations are decoupled.
The static problem splits up into (a) the tension/compression problem (the unknowns are: u, N , ε), (b) the bending in the plane x-z (the unknowns are: β, w, M y , T z , κ y , γ z , B z , γ z ), (c) the bending-torsion problem (the unknowns are: ϕ, v, θ, M z , T y , M, κ z , γ y , κ ω , ρ, B, B y , γ y ), (d) Simplifications due to bi-symmetry of the cross-section. Here, additionally, S 13 = 0, and the problem (c) above splits up into (d) the bending problem in the y-z plane (the unknowns are: ϕ, v, M z , T y , κ z , γ y , B y , γ y ), (e) the torsion problem (the unknowns are: θ, M, B, ρ, κ ω ). Thus, in the case considered, the problem of torsion is formulated in the same manner as in Vlasov's theory.

El Fatmi's general theory equipped with the warping functions of first-order approximation
The aim of the present section is to construct El Fatmi's [2,3] general theory (FG) by equipping it with the first-order warping functions (4.1).

The kinematic assumptions and the associated stress fields
The kinematic assumptions (5.1) are relaxed such that the non-uniform warping is modeled by independent variables ϑ, χ z , χ y : The primal unknowns of the theory form the collection u = u β w ϕ v θ ϑ χ z χ y . The fields ϑ, χ z , χ y are new primal unknowns, now unconstrained, i.e., not linked with the standard primal unknowns. The strains associated with the kinematic fields (6.1) are expressed by As before, the strains ε y , ε z , γ yz vanish, since within the plane y-z the cross-section A is subject only to translations and a rigid rotation. The nonzero stress components σ x , τ xy , τ xz are expressed by (1.2).

The equilibrium equations and the admissible boundary conditions
We start from computing the integrand of the virtual work of stresses associated with the strains (6.2): The expression above suggests introducing the stress resultants N , M y , M z , B, B y , B z , see (5.9), and new stress resultants: -the standard transverse forces -the axial moment of pure torsion Let us compute the virtual work of the loading. The work of the span load is, for simplicity, assumed as before (5.5). To fix attention, we assume that the end x = 0 is fully clamped, where all the components of the collection u = u , β , w , ϕ , v , θ , ϑ , χ z , χ y vanish. The end x = l is free and loaded by the tractions t x (y, z), t y (y, z), t z (y, z), see Sect. 5.2. We compute the virtual work of the tractions by (5.6), taking the virtual displacements in the form (6.1), to obtain L ends =Nū(l) +M yβ (l) +T zw (l) +Mθ(l) +Bθ(l) +B yχy (l) where the forces at the end x = l:N ,M y ,T z ,M,B,B y ,M z ,T y ,B z are tractions' resultants given by (5.8).
By equating the virtual work of internal forces (6.6) to the virtual work of the loads:L span +L ends , taking into account that the componentsū,β,w,φ,v,θ,θ,χ z ,χ y are arbitrary within (0,l] and vanish at x = 0, one obtains the equations of equilibrium similar to those found in Sect. 5 (they will not be displayed here), see (5.13), as well as the boundary conditions at the free end: Note that in this theory the effective end forces (like those in (5.15)) do not appear. Equation (5.4) determines the admissible boundary conditions of this theory. The conjugated pairs of the end-forces and the corresponding end-displacements are now: θ), (B, ϑ), (B z , χ z ), (B y , χ y ). (6.9) The number of standard admissible types of boundary conditions is 2 9 . For instance, the fork condition means prescribing   N , M y , w, M z , v, θ, B, B z , B y , (6.10) since then the tractions t x are applied, while the angle of torsion and the lateral displacements are fixed.

The constitutive equations
Substitution of the stress components into the definitions of the stress resultants (5.9), (6.5) and integration over the cross-section area leads to the constitutive equations of the form and C = C i j i, j=1,2,...,6 with the components given by C 44 = (ζ ,y ) 2 + (ζ ,z − 1) 2 , C 45 = ζ ,y (η ,y − 1) + η ,z (ζ ,z − 1) , C 46 = ζ ,y ω ,y + (ζ ,z − 1)ω ,z , the vector comprising the so-called primal internal transverse forces and the primal torsional moment: with the vector of strains due to shear and torsion: is diagonal, see the equation (37) in El Fatmi [2]. However, the choice of the warping functions assumed in the present paper (i.e., the first-order warping functions) does not lead to this diagonalization although it is the simplest acceptable choice. The general version of the El Fatmi theory loses its crucial property. The coupled form of the equations linking (6.16) with (6.17) makes the theory discussed here fairly complicated if used as an analytical tool.

Analysis of the static problem
The primal unknowns of the theory form the collection u = u β w ϕ v θ ϑ χ z χ y . Substitution of the constitutive equations into the virtual work equation (5.4) with the left-hand side given by (6.6) and the right-hand side given by (5.5) and (6.7) leads to the variational equation of a standard form: a F (u,ū) = f (ū) for eachū kinematically admissible (6.18) with a new bilinear form a F (u,ū) whose explicit definition is omitted here to save space. One can prove that this bilinear form is symmetric and positive definite, provided that the bar is correctly supported. This implies uniqueness of the solution u. To prove existence, more subtle regularity assumptions are necessary; these mathematical technicalities are omitted.

The theory of Kim and Kim (2005)
Kim and Kim [15] neglect the warping due to shear, since the theory is intended for modeling deformations of thin-walled bars. The convenient starting point to construct this theory is the general theory of El Fatmi, see Sect. 6. By assuming that the fields χ z , χ y vanish, and consequently, that χ z = 0, χ y = 0 one arrives at the Kim and Kim theory, originally proposed for the thin-walled profiles and extended here to bars of arbitrary cross-sections, possibly multi-connected. The primal unknowns of this theory form a column vector u = [u, β, w, ϕ, v, θ, ϑ] T . The theory is based on the kinematic assumptions (6.1) with neglecting the last two terms in the equation for the axial displacement u x . In Eqs. (6.2) defining the strains, we put: χ z = 0, χ y = 0, χ z = 0, χ y = 0. The quantities: (ε, κ z , κ y , ϑ , γ z , γ y , ρ, ϑ), where ϑ = dϑ/dx, play the role of strain measures of this theory, while the quantities (N , M z , M y , B, T z , T y , M ν M ω ) are the corresponding stress resultants. The set of the constitutive equations is composed of (6.11) and Unfortunately, the stiffnesses involved in (7.1) are overestimated, since the kinematic assumptions are chosen as too simple. Some stiffnesses should be made smaller by introducing shear corrections. To have better estimates of the stiffnesses in (7.1), Kim and Kim [15] make use of the saddle-point functional of Eric Reissner [20]; to construct the compliances of the theory, Kim and Kim apply the formulae for the stress distributions known from Vlasov's theory of thin-walled bars. Thus, the theory of Kim and Kim has been originally derived exclusively for thin-walled bars.
The left-hand side of the virtual work equation (5.4) takes now the form If the end x = 0 is fully clamped, there the displacement vector (7.1) vanishes, while at the free end x = l the components of this vector are arbitrary. The virtual work of the tractions at the free end equals L ends =Nū(l) +M yβ (l) +T zw (l) +M z (−φ(l)) +T yv (l) +M νθ (l) +Bθ(l) (7.3) where the end-forcesN , . . . ,B are given by (5.8). The virtual work of the span load is assumed as previously by (5.5). The equality (5.4) implies the variational equation of equilibrium (not displayed here to save space) as well as the boundary conditions at the free end: The equality (5.4) determines also the admissible boundary conditions within the theory. The conjugated pairs of end-forces and the corresponding end-displacements are listed below,

(N , u), (M y , β), (T z , w), (M z , −ϕ), (T y , v), (M ν , θ), (B, ϑ). (7.5)
Within this theory 2 7 types of standard boundary conditions are admissible. For instance, at the fork support the following quantities are prescribed : N , M y , w, M z , v, θ, B. We see that the conditions at the fork supports are the same as in Vlasov's theory.

The theory of Librescu and Song (2006)
Let us remind that the original Vlasov's theory applies to thin bars of thin-walled cross-sections in which transverse forces do not affect transverse shear. In Sect. 3 of the book by Librescu and Song [16], Vlasov's theory has been extended towards including the effects due to transverse shear deformation.

An attempt to construct the theory basing on Kim and Kim's theory
In the theory by Kim and Kim [15], the field ϑ is an independent variable. Within Vlasov's theory, this field is constrained by: ϑ = ρ = dθ/dx. This assumption leads to linking the torsional moments M ν and M ω , since J = C 33 + C 36 + C 63 + C 66 , see (5.18), (5.19). The constitutive equations (7.1) take the form ⎡ Thus, the above reduction to the Kim and Kim theory is unsuccessful, since the torsional moment, although referred to the shear/torsion center S, does depend upon the measures γ z , γ y as transverse shear. Moreover, upon inverting (8.2) we note that the measure of torsion does depend upon the transverse shear forces, although their lines of action cut the center S. We remember that within the El Fatmi simplified theory equipped with the same warping functions as used here the torsion has not been caused by the transverse forces, see (5.17). This test of correctness is not passed here. with the strain measures ( ε, κ y , γ z , κ z , γ y , ρ, κ ω ) according to (5.17): the first six equations along with the equation B = E I ω κ ω . The strains are expressed by the displacements u, β, w, ϕ, v, θ according to (5.2). Now the torsional moment M does not depend on γ z , γ y while the transverse forces T z , T y do not affect the measure of torsion ρ. The model takes into account correctly that the point S is the center of shear/torsion.

Construction based on the simplified version of El Fatmi's theory
The formal correctness of the theory is conditioned by the positive definiteness of the density energy W ; now the last term in (5.21) reduces to a single term E I ω (κ ω ) 2 /2. Similarly, in the definition of the bilinear form a F S (., .) one should simplify the last term of the integrand to the form E I ω d 2θ dx 2 d 2 θ dx 2 . The bilinear form (5.23) will preserve its positive properties, and hence the new version of the problem (5.24) is still well posed.

Vlasov-like theory
The point of departure will be the theory of Librescu and Song [16], Sect. 8. Vlasov's [1] theory concerns thin bars, and hence the transverse shear deformations are neglected. The Vlasov-like theory, applicable to both the solid and thin-walled profiles, will be based on puttingγ z = 0,γ y = 0 into (8.3). The virtual work of the end-loads at x = l and the virtual work of the span load are assumed as in Sect. 8. The primal unknowns are u, w, v, θ.
The variational equations of equilibrium take the form, the virtual work of the end forces being omitted for simplicity, 1.
where M eff = M − dB/dx is the effective torsional moment. The constitutive equations link the stress resultants N , M y M z , M, B with the strain measures ε,κ y ,κ z , ρ, κ ω defined by The transverse forces T z , T y are reactions of the constraints: β + dw dx = 0, ϕ + dv dx = 0. The density of energy (5.21) takes the form 5) and the bilinear form of the equilibrium problem reads now This bilinear form is symmetric and positive definite. The rigid body motions are the displacement fields for which the strain measures (9.3) vanish. They have the form: where c 1 , . . . , c 6 are constants. They vanish if the bar is correctly supported; then: a V (ū,ū) = 0 implies u = 0 , which assures uniqueness of the solution. The existence problem exceeds the framework of this paper.
The Vlasov-like theory involves only one stiffness: E I ω (the torsional warping stiffness) not present in Bernoulli-Euler's theory. Let us emphasize that in Vlasov's theory this stiffness is defined for thin-walled profiles. In the present paper, this stiffness is defined for all profiles, since the warping function ω (y, z) is determined by (2.2) by solving the elliptic problem (2.1). The torsional warping moment of inertia is given by I ω = (ω|ω); this formula is applicable to all shapes of A, as has been stressed already by Gruttmann et al. [10] and Wagner and Gruttmann [22].

The Timoshenko-like theory
The Timoshenko [23,24] theory of bars can be constructed by neglecting the warping due to torsion in the Librescu-Song theory. We putκ ω = 0 in Eq. The equilibrium problem splits up into the problems of: (a) tension/compression, (b) bending and transverse shear, (c) torsion. We note that: If the cross-section domain is asymmetric, then A yz = 0 and the bending in the plane x-z always induces bending in the plane y-z, and vice versa. This is the feature which distinguishes this theory from Vlasov's and Bernoulli-Euler's theories. Both the mentioned theories, as neglecting transverse shear effects, lead to a full decoupling between the bending phenomena in the principal planes.
The characteristic features of the Timoshenko-like theory are: there exist statically determinate problems in which the equilibrium equations can be solved first, paving the way to construct the diagrams of the strain measures and compute displacements by Maxwell-Mohr integration rule. Thus, the force method can be formulated and used successfully, which is crucial in the case of bars with stiffnesses varying along the bar axis, (f) in contrast to Bernoulli-Euler's theory, this theory includes length scales, e.g., E J y /G A y ; consequently, the numerical methods may suffer shear locking effects, (g) the primal unknowns can be (seemingly) interpreted in terms of the rigid body mechanics, since they represent the components of two vectors: of translation and rotation, (h) the stress resultants can also be (seemingly) interpreted in terms of the rigid body mechanics, since they represent the components of two vectors: of force and moment.
However, the stress resultants N , M y , M z and the primal variables u, β, ϕ are referred to the axis x linking the gravity centers, while the stress resultants T z , T y , M and the primal variables w, v, θ are referred to the axis x (s) linking the torsion/shear center, see Fig. 2. That is why both the interpretations (g), (h) are correct only if the cross-section is bisymmetric, since only then the gravity center and the shear/torsion center coincide. Thus, the commonly argued simplicity of interpretation of the primal variables and the stress resultants within Timoshenko's theory is vague.

The theory of Bernoulli-Euler
This theory can be constructed by putting I ω = 0 in Vlasov's equations. Then, the bimoment B vanishes, and the variational equation of equilibrium (9.1) is appropriately reduced by putting B = 0 in (9.14). The primal unknowns are u, w, v, θ. The conjugated pairs of end-forces and end-displacements are The constitutive equations linking the stress resultants N , M y M z , M with the strain measures ε,κ y ,κ z , ρ are fully decoupled if referred to the principal axes y, z, see Eqs. (9. 4.1-4). Thus, the static problem splits up into the problems of: (a) tension/compression, (b) bending in the plane x-z, (c) bending in the plane x-y, (d) torsion. The bending in the plane x-z never (i.e., independently of the shape of the cross-section) induces deflection in the y direction, and vice versa. Let us recall that Vlasov's theory has the same property, see Sect. 9. We conclude: Neglecting the transverse shear deformations implies decoupling between the bending phenomena in the two planes: x-z, y-z.
Just here we see a substantial difference between the theories of Vlasov and Bernoulli-Euler and the Timoshenko-like theory. This also arouses a suspicion that Vlasov's and Bernoulli-Euler's theories do not predict properly the specific behavior of bars of asymmetric cross-sections.
The main unknowns of all the theories discussed in the present paper are, for the convenience of the reader, set up in Table 1. The theories are also characterized by the conjugated pairs of the end-forces and the end-displacements, not displayed in Table 1 but given in the present paper in appropriate places. Table 1 shows that only in the theories by Timoshenko and Bernoulli-Euler the number of primal variables is equal to the number of stress resultants; hence only within these two theories there exists a class of statically determinate problems, and only within these theories the force method for static problems may be formulated.

Sectional characteristics of selected cross-sections
Application of the theories of bars is conditioned by having at our disposal the cross-sectional characteristics of the profiles. The crucial step is to construct the solutions ω 0 , ζ , η to the elliptic problems (2.1), (3.1), (3.3). All the computations have been performed with using the second author's program in C++ implementing the FEM algorithm of Section 4.4 and the program in Python which implements triangulation of the mesh, with using the methods: set_points(…), set_facets(…), build(…) of the package meshpy available in Python Environment in Visual Studio 2019. Moreover, the procedure SVD in C++ available in Press et al. [25] is used for solving the linear algebraic problem with a rectangular matrix. The warping functions are constructed for the case of Poisson's ratio ν = 0; the comparison of results found in other sources should take this assumption into account. For the selected ten profiles, we list here for future references all the characteristics involved in the theories discussed. The dimensions of all the cross-sections are given in (cm).
The coordinates (y 0 , z 0 ) of the gravity center O are computed in the initial Cartesian coordinate system y 0 , z 0 . The point O becomes the center of the Cartesian coordinate system y, z determined by the principal directions of the cross-section, and hence J yz = 0 in all cases; in the case of asymmetric cross-sections, the axis y is rotated by the angle α with respect to the axis cutting the point O and being parallel to y 0 , see Figs. 3, Table 1 Setting up the unknowns: the primal variables u, the strain measures ε, and the corresponding stress resultants σ of the discussed theories of straight elastic and prismatic bars El Fatmi's general theory 5,6,7,8,9,10,11, and 12. The given coordinates (y S , z S ) of the shear/torsion center S are referred to the y, z system. In Tables 2, 3, and 4, the dimensions of the sectional characteristics are all measured in cm n (n ≥ 2) or are non-dimensional

Librescu and Song's theory
The data and the results are given in the following sequence: the shape of the section A with dimensions in (cm), the sectional characteristics found in the literature or computed according to the formulae available in the literature. All the other necessary characteristics (computed by the present authors) are set up in Tables 3  and 4. Moreover, the E-notation: mE ± n = m 10 ±n is used to save space.

Example 1 Rectangular section (1R)
The shear correction factors of the rectangular section shown in Fig. 3, according to Table 2 in Gruttmann and Wagner [11], are k y = k z = 0.8333 (recall that we refer to the case of Poisson's ratio being zero). The torsion constant is J = 7.324 cm 4 , calculated analytically by the standard formula from the theory of pure torsion, see Sokolnikoff [26]. The necessary characteristics along with the components of the matrix S defined in (5.17) are given in Tables 3 and 4. Actually, the matrix S is diagonal.

Example 2 Ellipsoidal section (2E)
The elliptical section and its meshing are shown in Fig. 4. The well-known solution to problem (2.1) makes it possible to construct the exact formulae: for the profile in Fig. 4 the torsional characteristics are: J = 5.026548 cm 4 , I ω = 0.3769911 cm 6 . The shear correction factors according to Hutchinson [27], for ν = 0, are given by: Here k y = k ξ=1/2 = 0.807692, k z = k ξ=2 = 0.886364, and these analytical results practically coincide with our numerical results, see Table 3.  Actually, the matrix S is diagonal, see Table 4. It is worth noting that I ω for the rectangular cross-section is 3.51 times larger than I ω of the ellipsoidal cross-section inscribed into this rectangle.

Example 3 T-section (3S)
The shear correction factors for the cross-section of Fig. 5, according to Table 3 in Gruttmann and Wagner [11], are k y = 0.7395, k z = 0.6767. The components of the matrix S defined in (5.17) are given in Table 4. Actually, S 12 = S 23 = 0.

Example 5 Channel section (5C)
The shear correction factors of the profile in Fig. 7, according to Table 2 in Dong et al. [7], for ν = 0.3 read k y = 0.352, k z = 0.449, yet the thickness t may be a little different from that in [7], see Fig. 1, p. 1682 [7], where this data is missing.   Table 2 The rail 49E1 characteristics computed by SOFISTIK upon introducing the contour of the shape shown in Fig. 9 y 0 (cm) z 0 (cm) The profile of the rail 49E1 is approximated as in Fig. 9. Its characteristics are set up in Tables 2, 3, and 4. The components of the matrix S defined in (5.17) are given in Table 4. Actually, S 12 = S 23 = 0.

Example 8 Thin-walled I-section of unequal flanges (8I)
The components of the matrix S of the I-profile of Fig. 10 are given in Table 4. Actually S 12 = S 23 = 0. The torsional constant according to [1]: J = 1973.33 cm 4 . The warping moment of inertia, as computed by Vlasov's rule [1], is I ω = 2.47959 × 10 7 cm 6 , while according to Timoshenko-Gere [28,Appendix] is I ω = 2.57143 × 10 7 cm 6 , and the latter result is closer to that computed by our program, see Table 3.
Other necessary cross-sectional characteristics of the profiles of Figs 3,4,5,6,7,8,9,10,11, and 12 have been computed by the program developed by the second author; the results are set up in Tables 3 and 4.
The warping functions due to transverse shear involved in the kinematic assumptions (5.1), (6.1) correspond to the case of Poisson's ratio ν being zero. Due to this choice, the point S is simultaneously the center of shear and the center of torsion, which essentially simplifies the formulations of all the theories discussed. Now we shall study in brief the error introduced by neglecting Poisson's ratio effect while computing the shear correction factors. The analysis is confined to the two sections: the ellipse (2E, see Fig. 4) and the I-section of unequal flanges (8I, Fig. 10). Both the sections are symmetric with respect to the z axis, hence k yz = 0.  Table 4 The components of the matrix S of the profiles of Figs. 3, 4, 5, 6, 7, 8, 9, 10, 11, and 12 computed  In general, the shear correction factors do depend on ν. The algorithm of computing these coefficients (now denoted by k ν z , k ν y ) has been proposed in Ref. [11] and implemented by the second author of the present paper. As expected, Poisson's ratio affects to some extent the values of k ν y and less strongly the values of k ν z , characterizing the elliptical section, see Table 5. On the other hand, the Poisson ratio effect is negligible in the case of the thin-walled profile 8I, see the values set up in Table 6. Here the factors k ν z , k ν y are almost independent of ν. In both cases of the profiles, the Poisson ratio effect is visible but is not essential.

Selected illustrative problems of statics
The theories developed make it possible to analyze deformations of straight bars under static loads and describe the buckling (including the torsional and lateral buckling) phenomena by using the methods of Vlasov's theory of thin-walled bars, applicable now to bars of arbitrary cross-sections. Two problems below are selected: (a) the lateral bending due to vertical loads, reflecting consequences of asymmetry if the shear deformations are encompassed; here, the Timoshenko-like theory is used, and: (b) construction of the safety domain for the stability problem of a bar subject to pure bending in two orthogonal planes.
13.1 A cantilever subject to the transverse end-force, as predicted by the Timoshenko-like theory Let us consider a cantilever clamped at the end x = 0 and subject to a transverse force of magnitude P acting at the angle α to the axis y, applied at the shear/torsion center S, see Fig. 13a. The nonzero internal forces are: The nonzero strain measures have the form where the coefficients α zz /A, α zy /A = α yz /A, α yy A are the components of the matrix A −1 given by (4.17 Let us introduce the non-dimensional quantities The position of the point S upon deformation is given bȳ v S = b yy cos α + b yz sin α,w S = b zy cos α + b zz sin α (13.5) The points (v S ,w S ) lie on an ellipse given by (13.7) One of the principal axes of the ellipse makes an angle β to the axis y, and tan(2β) = 2b/(a − c) or It is seen that the angle β does not correspond to the principal direction of the matrix A, except for the rare case of J y = J z . If α yz = 0, then β = 0, and the principal axes of the ellipse coincide with the principal axes of the domain A. If the force P makes angles α = β, α = β + π/2, the magnitudes of displacements of the point S are extremal. The shape of the ellipse (13.6) reflects the geometric characteristics of the crosssection domain A and the slenderness of the bar. In general, the effect of a lateral deflection occurs to be marginal. For instance, if the bar is of channel-like cross section (5C) and of length l = 40 cm, in case of ν = 0.3, E = 210 GPa = 2.1 × 10 9 kg/(cm s 2 ) the angle of slope of the ellipse is small: β = 0.55442 • , hence not visible in Fig. 13b. The displacements are computed for P = 1.0 × 10 5 kg cm/s 2 .
13.2 Construction of the stability domain for the problem of a simultaneous action of two pairs of bending moments at the ends of the bar Consider a bar fork-supported at both the ends x = 0, x = l. The ends of the bars are subject to the external momentsM y ,M z directed opposite to the axis y and along the axis z, respectively, as shown in Fig. 14. This convention assures that the moments applied induce negative values of the stress σ x for positive z and y , respectively.
The subject of the discussion is the shape of the stability domain: F(M y ,M z ) ≤ 1. If the pair (M y ,M z ) satisfies the above condition, then the bar cannot suffer a lateral buckling. To define the function F, we introduce the auxiliary quantities: of length dimension. Let us introduce the referential forces  The stability domain is the interior of the circle given by see Vlasov [1] and Weiss and Giżejowski [29]. The position of this circle and its radius depend on the geometric characteristics of the cross-section A, on the material data, and on the length of the bar, see Fig. 14b. The results given in Tables 7, and 8 have been computed for the data: ν = 0.3, E = 210 GPa = 2.1 × 10 9 kg/(cm s 2 ), and for the subsequent values of the lengths l of the bar measured in (cm); the buckling load values (13.10) are measured in kg cm/s 2 , while other characteristics are non-dimensional.
Let us stress that the result (13.12) originally found by Vlasov within the theory of thin-walled profiles is here used for bars of sections of arbitrary shape, due to appropriate extension of the range of applicability of Vlasov's bar theory.

Final remarks
Prior to publication of the paper by El Fatmi [2,3], Vlasov's theory had been viewed as applicable only for thin-walled profiles. El Fatmi has shown that by appropriate construction of the kinematic assumptions one can formulate a theory which is a natural extension of Vlasov' s theory in two directions: to include warping due to transverse shear and to extend the theory to be applicable for bars of arbitrary shape of the cross-sections. The other result which has shed a new light on Vlasov's theory is due to Wagner and Gruttmann [22] who proved that the warping functions in Vlasov's theory are approximants of the warping functions of the theory of pure torsion. Just these ideas are the basis for the constructions of the bars' theories discussed in the present paper. The original Vlasov's theory is here treated as a mathematical scheme applicable for bars of arbitrary cross-sections, and also other theories of bars discussed are viewed as broadly applicable, including Kim and Kim's [15] and Librescu and Song's [16] approaches. The Vlasov-like theory discussed in Sect. 9 has the mathematical structure of the standard theory of thinwalled bars yet refers to bars of arbitrary cross-sections. The stiffnesses E I ω , E J are computed by (5.17)-(5.18) on the basis of the torsional warping function ω referred to the shear/torsion center. In the original Vlasov's theory, the warping moment of inertia I ω is computed as a contour integral along the middle axis of the bars' walls; here I ω is given by the integral over A, thus making this theory applicable to non-thin-walled profiles. This paves the way for extension of the renowned formulae known from Vlasov's theory for critical loads due to the lateral bucking and the lateral-twisting buckling towards predicting instable behavior of non-thin-walled bars of atypical shapes.
The 'raison d'être' of El Fatmi's general approach is the decoupling between the primal stress resultants (transverse forces and torsional moment) forming the vector (6.16) and the corresponding primal strain measures in (6.17). In the present paper, a proof is given that this decoupling does not hold if the warping functions are chosen as the first-order warping functions, see (4.1). Despite this controversy, the general El Fatmi theory is crucial as the natural starting point for constructing and justifying other theories of bars.
The assumed choice of the warping functions (4.1) leads to the coupled form of the constitutive equations corresponding to the transverse shear, see Eqs. (5.17.4,5). The effective areas (4.15) involved there are the same as those found by Schramm et al. [5] for the case of ν = 0. A proof is given that the relevant matrix (4.16) is symmetric and positive definite; the off-diagonal terms are nonzero if the cross-section is asymmetric, even if the axes y, z are principal axes of the cross-section. The algorithm of constructing this matrix is given with all details, including the numerical codes for solving the underlying scalar elliptic problems. This coupling phenomenon shows an essential difference between symmetric and asymmetric profiles and shows essential differences between the theories including transverse shear effects and those neglecting them. We conclude that neglecting the transverse shear deformations has far-reaching consequences of unjustifiedly neglecting the unavoidable lateral bending of asymmetric profiles. Moreover, taking into account three various cases (4.21), one obtains respectively. The process of computing the last (fourth) raw of the matrix K e is performed in a similar way as above and reduces to computing the three integrals for each finite element A e . All the integrals considered above have been computed analytically with using MAPLE system; the formulae defining the integrals have been rewritten to the form of a code in C (by calling the MAPLE's command C(expression, optimized, precision=double)) in order to implement the second author's program (in C++) of numerical construction of the solution to the Poisson equation (4.20.1) with the von Neumann boundary condition (4.20.2) augmented with the integral condition (4.20.3) and in order to compute all the cross-sectional integral characteristics involved in the discussed theories of bars.